Python-Continuous-time Markov chain state transition simulation

Overview

** Continuous Time Markov Chain (Continuous Time Markov Chain) state transition simulation is implemented. The algorithm used is the Gillespie Algorithm. After posting the implementation, I will briefly explain the outline.

** Discrete time ** For the Markov chain state transition simulation, refer to this article. Python --Markov chain state transition simulation

Implementation

Consider a Markov chain in the case of three states (state = 0,1,2).

CTMC.py

import numpy as np
import math
import random
import matplotlib.pyplot as plt
#Number of states
N=3
#Transition rate matrix
TRANSITION_MATRIX = [[-0.3, 0.1, 0.2],
                     [0.2, -0.9, 0.7],
                     [0.7, 0.2, -0.9]]

initial_state = 0 #initial state
current_state = initial_state #Current state
next_state = 0 #Next state
dt = 0 #State stay time

state_list = [] #State storage list
dt_list=[] #Storage list of staying time
state_list.append(current_state)
dt_list.append(dt)

time = 0
T=15 #Observation time

#Gillespie algorithm(Gillespie Algorithm)
while(time <= T):
    rand_s = random.random() #Random numbers that determine the state
    rand_t = random.random() #Random number that determines the length of stay
    q =  -TRANSITION_MATRIX[current_state][current_state] #propensity function
    total = 0
#Determine the next transition state
    for i in range(N):
        if(i == current_state):
            continue
        total += TRANSITION_MATRIX[current_state][i]
        if(rand_s < total/q):
            next_state  = i
            break;
    current_state = next_state
    state_list.append(current_state)
#Determine the length of stay in the current state
    dt = 1/q * math.log(1/rand_t)
    time += dt
    dt_list.append(dt)
#plot
arr_state = np.array(state_list)
arr_dt = np.array(dt_list)
arr_cumsum = np.cumsum(arr_dt)
arr_state_repeat = np.repeat(arr_state,2)
arr_time_repeat = np.repeat(arr_cumsum,2)
plt.plot(arr_time_repeat[1:],arr_state_repeat[:-1],label="state")
plt.xlabel("time")
plt.ylabel("state")
plt.title("DTMC")
plt.legend()
plt.grid(color='gray')

It is a graph of the result.

It can be seen that the transition is repeated between the states [0,1,2] until time T = 15. (In this graph, it transitions to around 17.5, but this is a specification. It can be removed by removing the end of the list, but it is not removed here because the code becomes complicated. )

Gillespie Algorithm

The main idea of this algorithm is according to the transition rate

1. When the next transition ** occurs **
2. ** Which ** state to transition to The two are randomly determined by random numbers. It is designed so that the larger the sum of the transition rates, the easier it is to transition immediately, and the larger the transition rate, the easier it is to transition. For more information, see Wikipedia.