How to implement a Markov Chain in python for a squared spatial grid?

Question

I'd like to understand how to implement in python a Markov Chain on a 2D grid (a 2D numpy array) given the values for all the states (positions) of the grid and also the transition matrix between all states.

Here I will use some figures to better explain my goal.

I have a 2x2 grid, each cell of the grid is identified by a state (the number of the grid cell) and by a value (some property of the cell).

I also created a possible transition matrix from each cell to each other cell of the grid. The value inside each cell of the transition matrix is the probability of transition from the cell selected by the column, to the cell selected by the row.

In this case for example the cell with state d can only transit to cell b, or c, with probability 0.6 and 0.4 respectively (that is Q(b,d) = 0.6 and Q(c,d) = 0.4).

At each transition I want to apply the function mean between the value stored in the cell at time t and the new values transiting in the cell at time t+1. For example after one cycle in this case (the transitions that actually took place in this case are presented under the arrow).

In each cell the old value has been averaged with the values of the cells that transitioned to that cell.

I'd like to create an efficient code that can implement this algorithm. I already have a code that does this but without markov chain and I think I may obtain significant better performance by using MC.

You have two different problems here, right? Computing the grid that results from a series of transitions doesn't require the Q matrix at all. — Tim Roberts, Jun 22 '23 at 17:40
you need the Q matrix to actually choose the transition to apply. Here I chose the transition manually (the ones below the arrow) but in a script you may want to do that with a rng. The transition `a->d` has been chosen because in the Q matrix `Q(d,a)!=0`. But you don't need the values of the actual probabilities in Q to explain the algorithm if that's what you meant, I added them so in case someone implement this on code already has a Q matrix to use. — Mirko, Jun 22 '23 at 17:43
"Markov chain" is just a model, not an algorithm. What you have described here IS a Markov chain. — Tim Roberts, Jun 22 '23 at 17:52
So are you asking "I have the transition matrix and the current state, how do I most efficiently compute a random next state using numpy?" — Stef, Jun 22 '23 at 17:55
This process with this Q matrix converges after about 40 generations. All the cells become equal. Is that what you expected? It surprised me. — Tim Roberts, Jun 22 '23 at 18:04
@TimRoberts i want to implement *an* algorithm that can do the things that i described, I do not know how to write code that does it efficiently. @ Stef Yes not necessarily just numpy, any other package is also fine. — Mirko, Jun 22 '23 at 18:06
Yes it should converge obviously, youre averaging all your cells at some point. Also consider that this happens because all cells are linked either directly or indirectly. In bigger grids with sparse transition matrixes (full of 0s) you may have many cells that are independent to one another and the averaging will highlight some cluster of connected cells. — Mirko, Jun 22 '23 at 18:06

score 0 · Answer 1 · answered Jun 22 '23 at 18:10

I'm going to throw this out here as a starting point. This uses numpy, but not for any matrix arithmetic, so improvements are possible. Note that I have implemented the original grid as a vector; there's no good reason to have it as a matrix, and you can always do print(value.reshape(2,2)) if you want to see it that way.

import numpy as np
import random

value = np.array([1,2,3,4]).astype(float)

Q = np.array([
    [ 0.0, 0.2, 0.3, 0.5 ],
    [ 0.5, 0.0, 0.0, 0.5 ],
    [ 0.1, 0.2, 0.0, 0.7 ],
    [ 0.0, 0.6, 0.4, 0.0 ]
])

# The transitions from the question.
# trans = ((0,3),(1,0),(2,3),(3,1))

def choose(Q):
    choices = []
    for fr,row in enumerate(Q):
        r = random.random()
        for to,val in enumerate(row):
            if r < val:
                choices.append( (fr,to) )
                break
            r -= val
    return choices

def update( old, trans ):
    new = old.copy()
    cnt = np.ones(value.shape,dtype=float)

    for fr,to in trans:
        new[to] += old[fr]
        cnt[to] += 1

    return new / cnt

for _ in range(100):
    value = update( value, choose(Q) )
    print(value)

How to implement a Markov Chain in python for a squared spatial grid?

1 Answers1