Hamzeh Alsalhi's blog

State Termination Probabilities in an Absorbing Markov Chain

Below is a succinct function solve that takes an absorbing Markov chain represented as a list of lists and returns a list of integers, one integer for each state, expressing the relative likelihood of terminating at that state.

Each state is represented as a row in the matrix. This row expresses the relative likelihoods of transitioning to each of the other states from this row. For example row [1, 0, 0, 0, 2] would mean we have a chance (33%) of transitioning to state 1 and a chance (66%) of transitioning to state 5.

To input a Markov chain edit the chain variable, keep in mind to stick to the conventions:

  • The first row is the start state
  • Have at least one terminal state
  • A terminal state can be written as a row of all zeros, or a row that has only a one in the column the same as the rows index. Example: A row five with all zeros and only a one in the fifth column would be a correct way to represent that state 5 in the absorbing Markov chain is a terminating state.

The solution is based around the idea that we can look at each pair of rows and reduce their function with regards to how they affect the probability of terminating at each state into one of the rows. There is no need for inverting or rearranging a matrix so it is an easy solution to code, the only function used that might not be considered a primitive in most languages is gcd for getting the greatest common divisor, but in most cases it is available and if not it is easily coded in a few extra lines.

import Data.List

chain = [[0, 2, 0, 3, 0], -- State 1) Starting Transient State
         [0, 0, 2, 2, 0], -- State 2) Transient State
         [0, 0, 1, 0, 0], -- State 3) Terminal State
         [0, 0, 0, 0, 1], -- State 4) Transient State
         [0, 0, 0, 0, 1]] -- State 5) Terminal State

solve :: [[Int]] -> [Int]
solve c =
    let m = [take i (c!!i) ++ [0] ++ drop (i + 1) (c!!i) | i <- [0..(length c) - 1]]
        terminal = filter (\r -> foldl (+) 0 (m!!r) == 0) [0..(length m) - 1]
        transient = [0..(length m) - 1] \\ terminal
    in foldl join m [(i,j) | j <- transient, i <- transient] !! 0

join m (r1, r2) =
    let sum_r2 = foldl (+) 0 (m!!r2)
        h1 i = (sum_r2 * m!!r1!!i) + (m!!r1!!r2 * m!!r2!!i)
        h2 i = if elem i [r1, r2] then 0 else h1 i
        new_row = map h2 [0..length (m) - 1]
        gd = foldl gcd (new_row!!0) new_row
        new_row_simplified = map (\x -> div x gd) new_row
    in take r1 m ++ [new_row_simplified] ++ drop (r1 + 1) m
    
main = print (solve chain)

Running this Haskell code gives

[0,0,1,0,4]

Which tells us there is no chance of terminating at state 1, 2, or 4 because 0 was returned in the position that maps to those states, this is correct because these are transient states. And there is a chance (20%) of terminating at state 3, and a chance (80%) of terminating at state 5.