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)

Question looks too general, ask Google about "coupling from the past".