Computing the Steady-State Vector of a Markov Chain - Maple Programming Help

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Applications and Example Worksheets : Statistics : examples/SteadyStateMarkovChain

Computing the Steady-State Vector of a Markov Chain

Introduction

This worksheet demonstrates the use of Maple to investigate Markov-chain models.

The following Maple techniques are highlighted:

• 

Creating a custom function

• 

Solving a specific example

• 

Creating a generic solution

Introduction to Markov Chains

A Markov Chain is a weighted digraph representing a discrete-time system that can be in any number of discrete states.  The nodes of the digraph represent the states, and the directed edge weight between two states a and b represents the probability (called the transition probability from a to b) that the system will move to state b in the next time period, given that it is currently in state a.  The sum of the transition probabilities out of any node is, by definition, 1.  The set of probabilities is stored in a transition matrix P, where entry (i, j) is the transition probability from state i to state j. Clearly, the sum of each row of P is 1.  

A well-known theorem of Markov chains states that the probability of the system being in state j after k time periods, given that the system begins in state i, is the (i, j) entry of Pk.

 

A common question arising in Markov-chain models is, what is the long-term probability that the system will be in each state?  The vector containing these long-term probabilities, denoted Pi, is called the steady-state vector of the Markov chain.

 

This Maple application creates a procedure for answering this question.  As a case study, we'll analyze a two-server computer network whose servers have known probabilities of going down or being fixed in any given hour. The goal is to compute the long-term probability that at least one server is working.

Algorithm for Computing the Steady-State Vector

We create a Maple procedure called steadyStateVector that takes as input the transition matrix of a Markov chain and returns the steady state vector, which contains the long-term probabilities of the system being in each state.  The input transition matrix may be in symbolic or numeric form.  

 

The procedure steadyStateVector implements the following algorithm:  Given an n x n transition matrix P, let I be the n x n identity matrix and Q = P - I.  Let e be the n-vector of all 1's, and b be the (n+1)-vector with a 1 in position n+1 and 0 elsewhere.  To compute the steady state vector, solve the following linear system for Pi, the steady-state vector of the Markov chain:

 

(Q | e)TPi=b

 

Appending e to Q, and a final 1 to the end of the zero-vector on the right-hand side ensures that the solution vector Pi has components summing to 1.

 

Procedure Code

Here is the steadyStateVector procedure.  The input is a transition matrix P, and the output is the steady-state vector pi reflecting the long-term probability of the system being in each state.  Comments are preceded by the # symbol.

 

steadyStateVector := proc( P::Matrix )

  #DECLARE LOCAL VARIABLES
  local n, Q, e, QT, b;
  
  #MAKE THE PROCEDURE SELF CONTAINED BY LOADING REQUIRED PACKAGES INSIDE THE PROCEDURE
  use LinearAlgebra in

  #EXTRACT THE DIMENSION OF THE TRANSITION MATRIX P
  n := Dimension( P )[1];

  #Q = P - I
  Q := P - IdentityMatrix( n );

  #e IS THE VECTOR OF ALL ONES
  e := <seq(1, i=1..n)>;

  #APPEND VECTOR e TO Q AND TRANSPOSE THE RESULT
  QT := Transpose( <Q | e> );

  #b IS THE UNIT VECTOR WITH 1 IN POSITION n+1
  b := UnitVector( n+1, n+1 );

  #SOLVES THE LINEAR SYSTEM QT*pi = b.
  return LeastSquares( QT, b );
  end use:
end proc:

Example Application: Reliability of a two-server network

The problem is to estimate the long-term probability that at least one server in a two-server computer network is working during any given hour.  We'll model this problem as a Markov chain as follows:  

 

Assume the network can be in one of three states:

1. 

Both servers are working

2. 

One server is working

3. 

Neither server is working  

 

Let:

• 

&lambda;1 be the probability that a server fails when both were okay an hour ago

• 

&lambda;2 be the probability that the second server fails when one was okay an hour ago

• 

&mu;1 be the probability that a broken server gets fixed when one was okay an hour ago

• 

&mu;2 be the probability that a broken server gets fixed when both were down an hour ago

 

The transition matrix is then the following:

 

P := Matrix( [ [1-mu[1],  mu[1],  0],
               [lambda[1],  1-(lambda[1]+mu[2]),  mu[2]],
               [0,  lambda[2],  1-lambda[2]] ] );

P1μ1μ10λ11λ1μ2μ20λ21λ2

 

Note that the steadyStateVector procedure computes pi symbolically.  Numerical values for lambda and mu are not required.

 

pi := steadyStateVector( P );

&pi;λ1λ2λ1λ2+λ2μ1+μ2μ1λ2μ1λ1λ2+λ2μ1+μ2μ1μ2μ1λ1λ2+λ2μ1+μ2μ1

 

pi is currently a symbolic vector. Let's supply some example values:

lambda[1]:=1/400; lambda[2]:=1/800; mu[1]:=1/4;  mu[2]:=1/8;

λ11400

λ21800

μ114

μ218

 

Ask Maple for the value of pi, and the updated vector is given, reflecting the inputs above.  By inspection, we see the vector pi is stochastic.

 

pi := steadyStateVector( P );

&pi;110101100101011000010101

 

The long-term probability that at least one server is operable in any given hour is the sum of the last two components of pi.

 

probWorking := pi[2] + pi[3];

probWorking1010010101

 

If we want a 10-digit floating-point approximation to this probability, use the evalf command.

 

evalf( probWorking, 10 );

0.9999009999

Erase the current values of lambda and mu.

lambda := 'lambda'; mu := 'mu';

λλ

μμ

 

Let's generalize this Markov chain, allowing for the possibility of both servers going down at once or both being repaired at once.

 

P := Matrix( [ [1-(mu[1]+mu[3]),  mu[1],  mu[3]],
               [lambda[1],  1-(lambda[1]+mu[2]),  mu[2]],
               [lambda[3],  lambda[2],  1-(lambda[2]+lambda[3])] ] );

P1μ1μ3μ1μ3λ11λ1μ2μ2λ3λ21λ2λ3

pi := steadyStateVector( P );

&pi;λ2λ1+λ3λ1+μ2λ3λ2λ1+λ3λ1+λ1μ3+μ1λ2+μ3λ2+μ1λ3+μ2λ3+μ2μ1+μ2μ3μ1λ2+μ3λ2+μ1λ3λ2λ1+λ3λ1+λ1μ3+μ1λ2+μ3λ2+μ1λ3+μ2λ3+μ2μ1+μ2μ3λ1μ3+μ2μ1+μ2μ3λ2λ1+λ3λ1+λ1μ3+μ1λ2+μ3λ2+μ1λ3+μ2λ3+μ2μ1+μ2μ3

 

What about a purely numeric matrix?

 

P2 := Matrix( [[1/4, 1/2, 1/4], [1/3, 1/3, 1/3], [1/4, 1/2, 1/4]]);

P2141214131313141214

pi := steadyStateVector( P2 );

&pi;273727

P3 := Matrix( [[.25, .45, .3], [.13, .33, .64], [.2, .6, .2]]);

P30.250.450.30.130.330.640.20.60.2

pi := steadyStateVector( P3 );

&pi;0.162983945355205880.44067176681185750.39569553441396604

 

Return to Index for Example Worksheets