Gillespie Algorithm在seir模型中的python实现

下面是使用 Gillespie Algorithm 实现 SEIR 模型的 Python 代码：

import numpy as np
import matplotlib.pyplot as plt


def gillespieSEIR(beta, sigma, gamma, N, E0, I0, R0, Tmax):
    # Initialize variables
    t = 0
    S, E, I, R = N - E0 - I0 - R0, E0, I0, R0
    S_list, E_list, I_list, R_list, t_list = [S], [E], [I], [R], [t]

    while t < Tmax:
        # Calculate rates of each event
        lambda1 = beta*S*I/N
        lambda2 = sigma*E
        lambda3 = gamma*I

        # Calculate total rate
        lambda_total = lambda1 + lambda2 + lambda3

        # Sample time until next event
        dt = np.random.exponential(scale=1/lambda_total)

        # Sample which event occurs
        event = np.random.choice(['lambda1', 'lambda2', 'lambda3'],
                                  p=[lambda1/lambda_total, lambda2/lambda_total, lambda3/lambda_total])

        # Update variables
        if event == 'lambda1':
            S -= 1
            I += 1
        elif event == 'lambda2':
            E -= 1
            I += 1
        else:
            I -= 1
            R += 1

        # Update time and lists
        t += dt
        S_list.append(S)
        E_list.append(E)
        I_list.append(I)
        R_list.append(R)
        t_list.append(t)

    return S_list, E_list, I_list, R_list, t_list


if __name__ == '__main__':
    # Set parameters
    beta = 0.2
    sigma = 0.1
    gamma = 0.05
    N = 1000
    E0, I0, R0 = 10, 1, 0
    Tmax = 100

    # Run simulation
    S, E, I, R, t = gillespieSEIR(beta, sigma, gamma, N, E0, I0, R0, Tmax)

    # Plot results
    plt.plot(t, S, label='Susceptible')
    plt.plot(t, E, label='Exposed')
    plt.plot(t, I, label='Infected')
    plt.plot(t, R, label='Recovered')
    plt.legend()
    plt.xlabel('Time')
    plt.ylabel('Number of individuals')
    plt.show()

其中，参数 `beta` 表示感染率，`sigma` 表示潜伏期转化为感染期的速率，`gamma` 表示康复率，`N` 表示总人口数，`E0`、`I0`、`R0` 表示初始时刻的潜伏期、感染期和康复人数，`Tmax` 表示模拟的时间长度。函数 `gillespieSEIR` 返回每个时间点的 S、E、I、R 数量，以及对应的时间点。最后，使用 matplotlib 库将结果可视化。