In a blog post I wrote in 2013, I showed how to simulate a discrete Markov chain. In this post we’ll (written with a bit of help from Geraint Palmer) show how to do the same with a continuous chain which can be used to speedily obtain steady state distributions for models of queueing processes for example.

A continuous Markov chain is defined by a transition rate matrix which shows the rates at which transitions from 1 state to an other occur. Here is an example of a continuous Markov chain:

This has transition rate matrix \(Q\) given by:

The diagonals have negative entries, which can be interpreted as a rate of no change. To obtain the steady state probabilities \(\pi\) for this chain we can solve the following matrix equation:

if we include the fact that the sum of \(\pi\) must be 1 (so that it is indeed a probability vector) we can obtain the probabilities in Sagemath using the following:

You can run this here (just click on ‘Evaluate’):

Q = matrix(QQ, [[-3, 2, 1], [1, -5, 4], [1, 8, -9]])
(transpose(Q).stack(vector([1,1,1])).solve_right(vector([0,0,0,1])))

This returns:

(1/4, 1/2, 1/4)

Thus, if we were to randomly observe this chain:

  • 25% of the time it would be in state 1;
  • 50% of the time it would be in state 2;
  • 25% of the time it would be in state 3.

Now, the markov chain in question means that if we’re in the first state the rate at which a change happens to go to the second state is 2 and the rate at which a change happens that goes to the third state is 1.

This is analagous to waiting at a bus stop at the first city. Buses to the second city arrive randomly 2 per hour, and buses to the third city arrive randomly 1 per hour. Everyone waiting for a bus catches the first one that arrives. So at steady state the population will be spread amongst the three cities according to \(\pi\).

Consider yourself at this bus stop. As all this is Markovian we do not care what time you arrived at the bus stop (memoryless property). You expect the bus to the second city to arrive 1/2 hours from now, with randomness, and the bus to the third city to arrive 1 hour from now, with randomness.

To simulate this we can sample two random numbers from the exponential distribution and find out which bus arrives first and ‘catch that bus’:

import random
[random.expovariate(2), random.expovariate(1)]

The above returned (for this particular instance):

[0.5003491524841699, 0.6107995795458322]

So here it’s going to take .5 hours for a bus to the second city to arrive, whereas it would take .61 hours for a bus to the third. So we would catch the bust to the second city after spending 0.5 hours at the first city.

We can use this to write a function that will take a transition rate matrix, simulate the transitions and keep track of the time spent in each state:

def sample_from_rate(rate):
    import random
    if rate == 0:
        return oo
    return random.expovariate(rate)

def simulate_cmc(Q, time, warm_up):
    Q = list(Q)  # In case a matrix is input
    state_space = range(len(Q))  # Index the state space
    time_spent = {s:0 for s in state_space}  # Set up a dictionary to keep track of time
    clock = 0  # Keep track of the clock
    current_state = 0  # First state
    while clock < time:
        # Sample the transitions
        sojourn_times = [sample_from_rate(rate) for rate in Q[current_state][:current_state]]
        sojourn_times += [oo]  # An infinite sojourn to the same state
        sojourn_times += [sample_from_rate(rate) for rate in Q[current_state][current_state + 1:]]

        # Identify the next state
        next_state = min(state_space, key=lambda x: sojourn_times[x])
        sojourn = sojourn_times[next_state]
        clock += sojourn
        if clock > warm_up:  # Keep track if past warm up time
            time_spent[current_state] += sojourn
        current_state = next_state  # Transition

    pi = [time_spent[state] / sum(time_spent.values()) for state in state_space]  # Calculate probabilities
    return pi

Here are the probabilities from the same Markov chain as above:

which gave (on one particular run):

[0.25447326473556037, 0.49567517998307603, 0.24985155528136352]

This approach was used by Geraint Palmer who is doing a PhD with Paul Harper and I. He used this to verify that calculations were being carried out correctly when he was trying to fit a model. James Campbell and I are going to try to use this to get an approximation for bigger chains that cannot be solved analytically in a reasonable amount of time. In essence the simulation of the Markov chain makes sure we spend time calculating probabilities in states that are common.