This is an archival dump of old wiki content --- see scipy.org for current material.

This page shows how to compute the stationary distribution pi of a large Markov chain. The example is a tandem of two M/M/1 queues. Generally the transition matrix P of the Markov chain is sparse, so that we can either use scipy.sparse or Pysparse. Here we demonstrate how to use both of these tools.

## Power Method

In this section we find pi by means of iterative methods called the Power method. More specifically, given a (stochastic) transition matrix P, and an initial vector pi_0, compute iteratively pi_n = pi_{n-1} P until the distance (in some norm) between pi_n and pi_{n-1} is small enough.

Fist we build the generator matrix Q for the related Markov chain. Then we turn Q into a transition matrix P by the method of uniformization, that is, we define P as I - Q/l, where I is the identity matrix (of the same size as Q) and l is the smallest element on the diagonal of Q. Once we have P, we approximate pi (the left eigenvector of P that satisfies pi = pi P) by the iterates pi_n = pi_0 P^n, for some initial vector pi_0.

More details of the above approach can be found in (more or less) any book on probability and Markov Chains. A fine example, with many nice examples and attention to the numerical solution of Markov chains, is `Queueing networks and Markov Chains' by G. Bolch et al., John Wiley, 2nd edition, 2006.

You can get the source code for this tutorial here: tandemqueue.py

```   1 #!/usr/bin/env python
2
3 labda, mu1, mu2 = 1., 1.01, 1.001
4 N1, N2 = 50, 50
5 size = N1*N2
6
7 from numpy import ones, zeros, empty
8 import scipy.sparse as sp
9 import  pysparse
10 from pylab import matshow, savefig
11 from scipy.linalg import norm
12 import time
```

This simple function converts the state (i,j), which represents that the first queue contains i jobs and the second queue j jobs, to a more suitable form to define a transition matrix.

```   1 def state(i,j):
2     return j*N1 + i
```

Build the off-diagonal elements of the generator matrix Q.

```   1 def fillOffDiagonal(Q):
2     # labda
3     for i in range(0,N1-1):
4         for j in range(0,N2):
5             Q[(state(i,j),state(i+1,j))]= labda
6     # mu2
7     for i in range(0,N1):
8         for j in range(1,N2):
9             Q[(state(i,j),state(i,j-1))]= mu2
10     # mu1
11     for i in range(1,N1):
12         for j in range(0,N2-1):
13             Q[(state(i,j),state(i-1,j+1))]= mu1
14     print "ready filling"
```

In this function we use scipy.sparse

```   1 def computePiMethod1():
2     e0 = time.time()
3     Q = sp.dok_matrix((size,size))
4     fillOffDiagonal(Q)
5     # Set the diagonal of Q such that the row sums are zero
6     Q.setdiag( -Q*ones(size) )
7     # Compute a suitable stochastic matrix by means of uniformization
8     l = min(Q.values())*1.001  # avoid periodicity, see the book of Bolch et al.
9     P = sp.speye(size, size) - Q/l
10     # compute Pi
11     P =  P.tocsr()
12     pi = zeros(size);  pi1 = zeros(size)
13     pi[0] = 1;
14     n = norm(pi - pi1,1); i = 0;
15     while n > 1e-3 and i < 1e5:
16         pi1 = pi*P
17         pi = pi1*P   # avoid copying pi1 to pi
18         n = norm(pi - pi1,1); i += 1
19     print "Method 1: ", time.time() - e0, i
20     return pi
```

Now use Pysparse.

```   1 def computePiMethod2():
2     e0 = time.time()
3     Q = pysparse.spmatrix.ll_mat(size,size)
4     fillOffDiagonal(Q)
5     # fill diagonal
6     x =  empty(size)
7     Q.matvec(ones(size),x)
8     Q.put(-x)
9     # uniformize
10     l = min(Q.values())*1.001
11     P = pysparse.spmatrix.ll_mat(size,size)
12     P.put(ones(size))
13     P.shift(-1./l, Q)
14     # Compute pi
15     P = P.to_csr()
16     pi = zeros(size);  pi1 = zeros(size)
17     pi[0] = 1;
18     n = norm(pi - pi1,1); i = 0;
19     while n > 1e-3 and i < 1e5:
20         P.matvec_transp(pi,pi1)
21         P.matvec_transp(pi1,pi)
22         n = norm(pi - pi1,1); i += 1
23     print "Method 2: ", time.time() - e0, i
24     return pi
```

Output the results.

```   1 def plotPi(pi):
2     pi = pi.reshape(N2,N1)
3     matshow(pi)
4     savefig("pi.eps")
5 if __name__ == "__main__":
6     pi = computePiMethod1()
7     pi = computePiMethod2()
8     plotPi(pi)
```

Here is the result:

## Improvements of this Tutorial

Include other methods such as Jacobi's method or Gauss Seidel.

SciPy: Cookbook/Solving_Large_Markov_Chains (last edited 2015-10-24 17:48:25 by anonymous)