ergodic markov chain - steady state

edit

A markov chain is just a stochastic matrix:

sage: m = matrix([[0.1, 0.2, 0.7], [0.5, 0.5, 0], [0.2, 0.2, 0.6]])
sage: m
[0.100000000000000 0.200000000000000 0.700000000000000]
[0.500000000000000 0.500000000000000 0.000000000000000]
[0.200000000000000 0.200000000000000 0.600000000000000]

As this one has positive entries, it is ergodic. Now, given a starting probability P0, you can iterate m on it:

sage: P0 = vector([0.2, 0.3, 0.5]) ; P0
(0.200000000000000, 0.300000000000000, 0.500000000000000)

sage: P = P0 * m^123 ; P               
(0.259740259740260, 0.285714285714286, 0.454545454545454)

Now, if you want to sample according to this probability vector, you can make 3 boxes of the interval [0,1] whose length correspond to your probability vector, pick a uniformly random point in [0,1], and check in which box it falls:

sage: def sample(P):
....:     s = 0
....:     a = random()
....:     for box in range(len(P)):
....:         s += P[box]
....:         if a <= s:
....:             return box
sage: sample(P)
1
sage: sample(P)
2
sage: sample(P)
0
sage: sample(P)
0
sage: sample(P)
1
sage: sample(P)
2

Looking for the steady-state probability, you can approach it by iterating the matrix m a lot:

sage: vector([1,0,0]) * m^1000
(0.259740259740259, 0.285714285714285, 0.454545454545454)

The convergence (quality and speed) is guaranteed by Perron Frobenius theorem.

If you are not sure how big should be the exponend to reach the steady-state, you can let Sage discover a fixed point by itself:

sage: def steady_state_iteration(m):
....:     # just in case vector like (1,0,0,0) falls in the contracting space.
....:     a = random()
....:     V = vector([a, 1-a] + [0 for i in range(len(m[0])-2)])
....:     while True:
....:         W = V * m
....:         W = W/sum(W)
....:         if V == W:
....:             return V
....:         else:
....:             V = W

sage: steady_state_iteration(m)
(0.25974025974, 0.285714285714, 0.454545454545)

You can also ask Sage for the dominating eigenvector. But you will notice that Sage does not accept to compute eigenvectors on the Real Field because no efficient and reliable was written there:

sage: m.left_eigenvectors()
NotImplementedError: eigenspaces cannot be computed reliably for inexact rings such as Real Field with 53 bits of precision,
consult numerical or symbolic matrix classes for other options

A first possibility is to work on an exact field like the rationals, and look for the (normalized) eigenvector whose eigenvalue is 1:

sage: def steady_state_eigenvector(m):
....:     for eig in m.left_eigenvectors():
....:         if eig[0] == 1:
....:             v = eig[1][0]
....:             return v/sum(v)

sage: m = matrix([[1/10,2/10,7/10],[5/10,5/10,0],[2/10,2/10,6/10]]) ; m
[1/10  1/5 7/10]
[ 1/2  1/2    0]
[ 1/5  1/5  3/5]

sage: steady_state_eigenvector(m)
(20/77, 2/7, 5/11)
sage: steady_state_eigenvector(m).change_ring(RR)
(0.259740259740260, 0.285714285714286, 0.454545454545455)

So, you can check it is quite close to the numerical approach above.

Another possiblity (suggested by @ppurka in comments), is to work in the ... (more)