Monads for Markov chains

Suppose you need to model a finite Markov chain in code. There are essentially two ways of doing that: one is to simply run a simulation of the Markov chain using a random number generator to obtain dice rolls and random cards from the decks, the other is to create a stochastic matrix containing the transition probabilities for each pair of states. In this post I will show how a single monadic description of the Markov chain dynamics can be used to obtain both a simulator and the transition matrix.

Let’s start with an example of Markov chain and how we would like to be able to implement in Haskell. Consider a simplified version of the familiar Monopoly game: there are 40 squares (numbered 0 to 39), you throw two 6-sided dice each turn, some special squares have particular effects (see below), if you get a double roll three times in a row, you go to jail. The special squares are: 30: go to jail 2, 17, 33: Community Chest 7, 22, 36: Chance Community Chest (CC) and Chance (CH) make you take a card from a deck and move to some other place depending on what’s written on the card. You will find the details on the code, so I won’t explain them here. This is of course a Markov chain, where the states can be represented by:

and a description of the game can be given in a monadic style like this:

As you can see, Game is a state monad, with an additional function choose that gives us a random element of a range:

This can be implemented very easily using the (strict) state monad and a random generator:

The runSim function runs the simulation and returns the list of visited states. This is already quite nice, but the best thing is that the same code can be used to create the transition matrix, just swapping in a new implementation of the Game type alias:

The idea is that we keep track of all possible destination states for a given state, with associated conditional probabilities. For those familiar with Eric Kidd’s series on probability monads, this is basically:

Now, how to get a transition matrix from such a monad? Of course, we have to require that the states are indexable:

So we iterate over all states and use the probability values contained in the monad to fill in the array cells corresponding to the selected state pair. To actually apply this to our Monopoly example, we need to make GameState indexable:

then finally we can try:

where .* is a simple vector-matrix multiplication function:

Calling statDistr 100 will return an association list of states with corresponding probability in an approximation of the stationary distribution, computed by applying the power method to the transition matrix. The number 100 is a pure guess, I don’t know how to estimate the number of iterations necessary for convergence, but that is out of the scope of this post, anyway.