Simulation de transition d'état de chaîne de Markov en temps continu Python

Aperçu

** Temps continu ** Implémente une simulation de transition d'état de la chaîne de Markov. L'algorithme utilisé est l'algorithme de Gillespie. Après avoir publié la mise en œuvre, j'expliquerai brièvement le plan.

** Temps discret ** Pour la simulation de transition d'état de la chaîne de Markov, veuillez vous référer à cet article. Python - Simulation de transition d'état de la chaîne de Markov

la mise en oeuvre

Considérons la chaîne de Markov dans le cas de trois états (état = 0,1,2).

CTMC.py

import numpy as np
import math
import random
import matplotlib.pyplot as plt
#Nombre d'états
N=3
#Matrice des taux de transition
TRANSITION_MATRIX = [[-0.3, 0.1, 0.2],
                     [0.2, -0.9, 0.7],
                     [0.7, 0.2, -0.9]]

initial_state = 0 #Etat initial
current_state = initial_state #État actuel
next_state = 0 #État suivant
dt = 0 #État du temps de séjour

state_list = [] #Liste de stockage d'état
dt_list=[] #Liste de stockage du temps de séjour
state_list.append(current_state)
dt_list.append(dt)

time = 0
T=15 #Temps d'observation

#Algorithme de Gillespy(Gillespie Algorithm)
while(time <= T):
    rand_s = random.random() #Nombre aléatoire qui détermine l'état
    rand_t = random.random() #Nombre aléatoire qui détermine la durée du séjour
    q =  -TRANSITION_MATRIX[current_state][current_state] #propensity function
    total = 0
#Déterminer le prochain état de transition
    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)
#Déterminer la durée du séjour dans l'état actuel
    dt = 1/q * math.log(1/rand_t)
    time += dt
    dt_list.append(dt)
#terrain
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')

C'est un graphique du résultat.

On voit que la transition se répète entre les états [0,1,2] jusqu'au temps T = 15. (Dans ce graphique, il passe à environ 17,5, mais c'est une spécification. Il peut être supprimé en supprimant la fin de la liste, mais il n'est pas supprimé ici car le code devient compliqué. )

Algorithme de Gillespie

L'idée principale de cet algorithme est en fonction du taux de transition

1. Quand la prochaine transition ** se produit **
2. ** Vers quel ** état effectuer la transition Les deux sont déterminés au hasard par des nombres aléatoires. Il est conçu de telle sorte que plus la somme des taux de transition est élevée, plus il est facile de faire la transition immédiatement, et plus le taux de transition est élevé, plus la transition est facile. Pour plus d'informations, veuillez consulter Wikipedia.