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

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 ${P}^{k}$.   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 $\mathrm{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 $\mathrm{Pi}$, the steady-state vector of the Markov chain:   ${\mathrm{\left(Q | e\right)}}^{T}\mathrm{Pi}=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 $\mathrm{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 $\mathrm{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 := ;   #APPEND VECTOR e TO Q AND TRANSPOSE THE RESULT   QT := Transpose( );   #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:

 • ${\mathrm{λ}}_{1}$ be the probability that a server fails when both were okay an hour ago
 • ${\mathrm{λ}}_{2}$ be the probability that the second server fails when one was okay an hour ago
 • ${\mathrm{μ}}_{1}$ be the probability that a broken server gets fixed when one was okay an hour ago
 • ${\mathrm{μ}}_{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]] ] );

${P}{≔}\left[\begin{array}{ccc}{1}{-}{{\mathrm{\mu }}}_{{1}}& {{\mathrm{\mu }}}_{{1}}& {0}\\ {{\mathrm{\lambda }}}_{{1}}& {1}{-}{{\mathrm{\lambda }}}_{{1}}{-}{{\mathrm{\mu }}}_{{2}}& {{\mathrm{\mu }}}_{{2}}\\ {0}& {{\mathrm{\lambda }}}_{{2}}& {1}{-}{{\mathrm{\lambda }}}_{{2}}\end{array}\right]$

Note that the steadyStateVector procedure computes $\mathrm{pi}$ symbolically.  Numerical values for $\mathrm{lambda}$ and $\mathrm{mu}$ are not required.

 > pi := steadyStateVector( P );

${\mathrm{π}}{≔}\left[\begin{array}{c}\frac{{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}}{{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}}\\ \frac{{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}}{{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}}\\ \frac{{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}}{{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}}\end{array}\right]$

$\mathrm{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;

${{\mathrm{\lambda }}}_{{1}}{≔}\frac{{1}}{{400}}$

${{\mathrm{\lambda }}}_{{2}}{≔}\frac{{1}}{{800}}$

${{\mathrm{\mu }}}_{{1}}{≔}\frac{{1}}{{4}}$

${{\mathrm{\mu }}}_{{2}}{≔}\frac{{1}}{{8}}$

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

 > pi := steadyStateVector( P );

${\mathrm{π}}{≔}\left[\begin{array}{c}\frac{{1}}{{10101}}\\ \frac{{100}}{{10101}}\\ \frac{{10000}}{{10101}}\end{array}\right]$

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

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

${\mathrm{probWorking}}{≔}\frac{{10100}}{{10101}}$

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';

${\mathrm{\lambda }}{≔}{\mathrm{\lambda }}$

${\mathrm{\mu }}{≔}{\mathrm{\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])] ] );

${P}{≔}\left[\begin{array}{ccc}{1}{-}{{\mathrm{\mu }}}_{{1}}{-}{{\mathrm{\mu }}}_{{3}}& {{\mathrm{\mu }}}_{{1}}& {{\mathrm{\mu }}}_{{3}}\\ {{\mathrm{\lambda }}}_{{1}}& {1}{-}{{\mathrm{\lambda }}}_{{1}}{-}{{\mathrm{\mu }}}_{{2}}& {{\mathrm{\mu }}}_{{2}}\\ {{\mathrm{\lambda }}}_{{3}}& {{\mathrm{\lambda }}}_{{2}}& {1}{-}{{\mathrm{\lambda }}}_{{2}}{-}{{\mathrm{\lambda }}}_{{3}}\end{array}\right]$

 > pi := steadyStateVector( P );

${\mathrm{π}}{≔}\left[\begin{array}{c}\frac{{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{3}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\lambda }}}_{{3}}}{{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{3}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\mu }}}_{{3}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{3}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{3}}}\\ \frac{{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{3}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{3}}}{{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{3}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\mu }}}_{{3}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{3}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{3}}}\\ \frac{{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\mu }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{3}}}{{{\mathrm{\lambda }}}_{{2}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{3}}{}{{\mathrm{\lambda }}}_{{1}}{+}{{\mathrm{\lambda }}}_{{1}}{}{{\mathrm{\mu }}}_{{3}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{3}}{}{{\mathrm{\lambda }}}_{{2}}{+}{{\mathrm{\mu }}}_{{1}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\lambda }}}_{{3}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{1}}{+}{{\mathrm{\mu }}}_{{2}}{}{{\mathrm{\mu }}}_{{3}}}\end{array}\right]$

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]]);

${\mathrm{P2}}{≔}\left[\begin{array}{ccc}\frac{{1}}{{4}}& \frac{{1}}{{2}}& \frac{{1}}{{4}}\\ \frac{{1}}{{3}}& \frac{{1}}{{3}}& \frac{{1}}{{3}}\\ \frac{{1}}{{4}}& \frac{{1}}{{2}}& \frac{{1}}{{4}}\end{array}\right]$

 > pi := steadyStateVector( P2 );

${\mathrm{π}}{≔}\left[\begin{array}{c}\frac{{2}}{{7}}\\ \frac{{3}}{{7}}\\ \frac{{2}}{{7}}\end{array}\right]$

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

${\mathrm{P3}}{≔}\left[\begin{array}{ccc}{0.25}& {0.45}& {0.3}\\ {0.13}& {0.33}& {0.64}\\ {0.2}& {0.6}& {0.2}\end{array}\right]$

 > pi := steadyStateVector( P3 );

${\mathrm{π}}{≔}\left[\begin{array}{c}{0.16298394535520588}\\ {0.4406717668118575}\\ {0.39569553441396604}\end{array}\right]$

 >