将 ODE 模型拟合到基于代理的模拟的感染数据

Question

我正在尝试使用 scipy.optimize 将 Susceptible-Infectious-Recovered (SIR) ODE 模型拟合到 python 中的感染数据，但是 ODE 模型的感染曲线导致延迟上升感染和高估的峰高（参见下面的橙色线）。

数据来自基于代理的模拟，具有均匀混合和三种状态（S、I 和 R）。理论上，这应该符合ODE模型的假设，并且模型应该能够完美拟合。我不确定我的拟合方法是否需要改进，或者基于代理的模拟是否创建了无法使用简单 ODE 模型重现的形状。

以下是数据和我的尝试。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

inf_data = [0.005,0.007,0.009,0.012,0.015,0.019,0.024,0.03,0.036,0.044,0.053,0.063,0.074,0.086,0.098,0.112,0.126,0.141,0.157,0.172,0.188,0.203,0.218,0.232,0.247,0.259,0.272,0.282,0.292,0.302,0.308,0.315,0.318,0.321,0.322,0.321,0.319,0.317,0.312,0.307,0.301,0.294,0.286,0.279,0.27,0.262,0.253,0.244,0.236,0.226,0.217,0.207,0.197,0.188,0.18,0.171,0.163,0.155,0.147,0.14,0.132,0.125,0.119,0.112,0.106,0.1,0.095,0.089,0.085,0.08,0.076,0.071,0.067,0.064,0.06,0.056,0.053,0.05,0.047,0.044,0.041,0.039,0.037,0.034,0.032,0.03,0.028,0.026,0.025,0.023,0.022,0.021,0.019,0.018,0.017,0.016,0.015,0.014,0.013,0.012,0.011] #percentage of infected
time_points = np.arange(len(inf_data))

def sir_model(y, x, beta, gamma):
    S, I, R = y
    dS = -beta * S * I / N
    dI = beta * S * I / N - gamma * I
    dR = gamma * I
    return dS, dI, dR

def fit_odeint(time_points, beta, gamma):
    return integrate.odeint(sir_model, (S0, I0, R0), time_points, args=(beta, gamma))[:,1]

# Initial values
N = 1.0
I0 = inf_data[0]
S0 = N - I0
R0 = 0.0

popt, pcov = optimize.curve_fit(fit_odeint, time_points, inf_data)
fitted = fit_odeint(time_points, *popt)

plt.plot(time_points, inf_data, 'o')
plt.plot(time_points, fitted)
plt.show()

Answer 1

我认为你的模型方程没有足够的“发挥”来拟合这些数据。作为曲线拟合的练习（如果不是在感染行为中），您可以在 dI/dt 的源项中添加另一个参数 p（以及 dS/dt 的下沉）。

请注意，R 在您的方程组中是多余的（因为 S+I+R=N） - 所以我在下面删除了它。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

inf_data = [0.005,0.007,0.009,0.012,0.015,0.019,0.024,0.03,0.036,0.044,0.053,0.063,0.074,0.086,0.098,0.112,0.126,0.141,0.157,0.172,0.188,0.203,0.218,0.232,0.247,0.259,0.272,0.282,0.292,0.302,0.308,0.315,0.318,0.321,0.322,0.321,0.319,0.317,0.312,0.307,0.301,0.294,0.286,0.279,0.27,0.262,0.253,0.244,0.236,0.226,0.217,0.207,0.197,0.188,0.18,0.171,0.163,0.155,0.147,0.14,0.132,0.125,0.119,0.112,0.106,0.1,0.095,0.089,0.085,0.08,0.076,0.071,0.067,0.064,0.06,0.056,0.053,0.05,0.047,0.044,0.041,0.039,0.037,0.034,0.032]
time_points = np.arange(len(inf_data))

def sir_model(y, x, beta, gamma, p):
    S, I = y
    dS = -beta * S * I ** p / N
    dI = beta * S * I ** p / N - gamma * I
    return dS, dI

def fit_odeint(time_points, beta, gamma,p):
    return integrate.odeint(sir_model, (S0, I0), time_points, args=(beta, gamma,p))[:,1]

# Initial values
N = 1.0
I0 = inf_data[0]
S0 = N - I0

popt, pcov = optimize.curve_fit(fit_odeint, time_points, inf_data)
fitted = fit_odeint(time_points, *popt)
print(*popt)

plt.plot(time_points, inf_data, 'o')
plt.plot(time_points, fitted)
plt.show()

您可能需要使用边界来删除警告，但它在我的机器上运行以给出

beta, gamma, p = 0.17905033272107485 0.06899938024967406 0.8725116899202809

请注意，幂 p 略小于 1：您需要初始感染率上升得更快一些。