Find probability of ending up in a row of a stochastic matrix

markov chainsmatricesprobabilitystochastic-matrices

I'm solving a programming problem which (after a small modification) turned out that can be represented as a stochastic matrix. One example is:

$$
\begin{bmatrix}
0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\
\frac{4}{9} & 0 & 0 & \frac{3}{9} & \frac{2}{9} & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}
$$

You can see there are $4$ "terminal states", represented by the 3, 4, 5 and 6 rows. As in, there's no way out from them.

The question then is "starting at row 1, calculate probability of the terminal states" and the answer to this particular matrix is:

$$
\begin{bmatrix}
0 \\
\frac{3}{14} \\
\frac{2}{14} \\
\frac{9}{14}
\end{bmatrix}
$$

Before I realized this is a stochastic matrix, I tried to implement some crazy complicated algorithm, where I would create a tree of all possible transitions, then I would go along every possible path and try to create a formula for each one; the tree also has cycles that need to be treated specially in the formula (you get a recursive probability), all of that is pretty complicated.

But now that I know it's a common concept called "stochastic matrix", I'm wondering, isn't there some common algorithm for calculating the probabilities of getting into the terminal states?

I can't seem to find anything related. It's either "online calculators" (which doesn't even give the probability of ending up in a row), or calculations of probability of being in a particular row at n-th step, which is something I'm not interested in.

Best Answer

@user8675309 gave helpful hints in the comments: on the naming of the problem and also that more explanation on these concepts may be found in a free book "Grinstead and Snell’s Introduction to Probability", Chapter 11. Most of the stuff that follows I obtained from the book.

The kind of stochastic matrices with a terminal state (a row) is called Absorbing Markov Chains, and the term for a "terminal state" is absorbing state. A non-absorbing state is called transient state. I'll also remark that every row in a stochastic matrix is a "state" from an according Markov chain, and the matrix represents probability of moving from one state to another. So my question may also be re-worded as "Probability of ending up in an absorbing state of a stochastic matrix". (not gonna change the title though, because that's currently a query a people unfamiliar with these concepts may ask)

Using this new terminology we can see that the matrix from the question has 2 transient states:

$$ \begin{bmatrix} 0 & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\ \frac{4}{9} & 0 & 0 & \frac{3}{9} & \frac{2}{9} & 0 \end{bmatrix} $$

and 4 absorbing states

$$ \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} $$

Then, steps to solving this problem are:

  1. Canonicalize the matrix. This means: reorder rows so the transient states go first and absorbing ones next. (strictly speaking, for the purpose of our problem you can have them any order you want as long as transient states are one by one. That is because in our calculation only the transient states will be used. To confirm that I tested different configurations of the matrix in question and they all gave me expected results)

    While moving rows around you have to rewrite references to them from other rows, but the algorithm is actually simple:

    1. Find two rows you want to swap and remember their indices. Let's say it's $i$ and $j$
    2. Swap the rows
    3. For each row in the matrix, swap i-th and j-th cells. So basically, you end up swapping i-th and j-th columns. Example code in python: 
      def swap(obj, i, j):
    tmp = obj[i]
    obj[i] = obj[j]
    obj[j] = tmp

def swap_rows(m, i, j):
    swap(m, i, j)
    # now fix references to the swapped rows
    for i_row, row in enumerate(m):
        swap(row, i, j)

    The matrix in the question is already in a canonical form, so no additional work is needed.

  2. Break the transient rows into matrices $R$ and $Q$ as follows: note that the matrix has an identity submatrix in the bottom-right corner. Let's call its size $t×t$. Then, $R$ will be the upper-right submatrix of transient states of size $r×t$, where $r$ are the rows from the top that are not included into the $t×t$. Then the $Q$ is the $r×r$ submatrix on the left of the $R$.

    In the question the submatrices are:

    $$ Q = \begin{bmatrix} 0 & \frac{1}{2} \\ \frac{4}{9} & 0 \end{bmatrix} R = \begin{bmatrix} 0 & 0 & 0 & \frac{1}{2} \\ 0 & \frac{3}{9} & \frac{2}{9} & 0 \end{bmatrix} $$

  3. Find an identity matrix $I_Q$ of $Q$

  4. Calculate "fundamental matrix" $N = (I_Q - Q)⁻¹$

  5. Multiply $NR$ matrices. The first row of the result will be the answer.

This will give the expected $0, \frac{3}{14}, \frac{1}{7}, \frac{9}{14}$.

In terms of a bonus, here's a dirty code to calculate the matrix in question with python numpy. Note, that numpy unfortunately doesn't work well with fractions module at the moment, so this example calculates result as a floating point (which sums up to 1, barring rounding errors).

import numpy as np

def find_fst_absorbing_row(m):
    for i_row, row in enumerate(m):
        if row[i_row] == 1 and sum(row) == 1:
            return i_row

m = [[0,   1/2, 0, 0,   0,   1/2],
     [4/9, 0,   0, 3/9, 2/9, 0],
     [0,   0,   1, 0,   0,   0],
     [0,   0,   0, 1,   0,   0],
     [0,   0,   0, 0,   1,   0],
     [0,   0,   0, 0,   0,   1]]

fst_absorbing_row = find_fst_absorbing_row(m)

Q = []
for i in range(0, fst_absorbing_row):
    Q.append(m[i][:fst_absorbing_row])
R = []
for i in range(0, fst_absorbing_row):
    R.append(m[i][fst_absorbing_row:])

Iq = np.identity(len(Q))
N = np.linalg.inv(Iq - Q)
result = np.dot(N, R)[0]

print(f'Expected: [0 0.21428571 0.14285714 0.64285714]\nActual: {result}')

I have a code that calculates all that stuff in fractions. For purposes of my programming challenge I had to re-implement matrix arithmetic without any libraries, and fractions module worked this way without any special modifications at all. So, basically, if you want fractions, you may use the standard fractions python module but without numpy.

Possible gotchas during implementation:

  • make sure a single-element matrix is handled [[0]]
  • if you had to canonicalize the matrix, make sure that upon returning the result the probabilities in it are sorted the way states were in the original (i.e. before canonicalization happened) matrix.
Related Question