SIR模型曲线拟合效果不佳问题排查与新数据预测咨询

拟合效果差的原因 #

1. 初始值设置错误 #

你的代码里重复定义了初始值：先设置I0=0.3，后面又覆盖为I0=ydata[0]（也就是0.0）。经典SIR模型中，若初始感染人数I0=0，则没有传染源，模型无法启动传播过程，拟合自然失效。

2. 数据不符合SIR模型假设 #

从数据和拟合图来看，感染比例出现多次起伏（比如第18天达到1.0后又回落，后续再次上升），但经典SIR模型的核心假设是：人群分为易感者(S)、感染者(I)、康复者(R)，康复者获得永久免疫，不会再次感染。因此感染人数应该呈现「单峰上升后单调下降」的趋势，你的数据波动明显违背了这个假设，可能的原因包括：

• 数据是新增感染数而非累计感染比例，但你用SIR模型拟合的是累计感染的比例；
• 数据存在噪声或统计误差；
• 实际场景中存在重复感染、防控措施变化（如封控、放开）等，经典SIR模型无法覆盖这类复杂情况。

3. 模型选择不匹配 #

如果实际场景存在重复感染，应该使用SIRS模型（康复后失去免疫，重新回到易感人群）；如果有潜伏期，考虑SEIR模型；如果有干预措施，需要给传播率beta加入时间依赖项（比如分段函数）。

修正后的拟合代码 #

先修正初始值问题，同时调整初始感染人数为非零值，示例如下：

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

y_total = [0.0, 0.0010131712259371835, 0.0035460992907801418, 0.00911854103343465,
    0.008611955420466059, 0.021783181357649443, 0.00911854103343465, 0.07852077001013172, 0.4397163120567376, 
    0.21681864235055726, 0.232016210739615, 0.5278622087132725, 0.13576494427558258, 0.2988855116514691, 0.37436676798378926, 
    0.4209726443768997, 0.544579533941236, 0.7254305977710233, 1.0, 0.7740628166160081, 0.43617021276595747, 0.48226950354609927]

x_total = range(0,22)
ydata = np.array(y_total, dtype=float)
xdata = np.array(x_total, dtype=float)

# 修正初始值：取第一个非零的感染值作为I0，避免初始值为0
I0 = ydata[ydata > 0][0] if any(ydata > 0) else 1e-6
N = 1.0
S0 = N - I0
R0 = 0.0

def sir_model(y, x, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

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

# 给curve_fit设置初始猜测值，提升拟合稳定性
popt, pcov = optimize.curve_fit(fit_odeint, xdata, ydata, p0=[1.0, 0.1])
fitted = fit_odeint(xdata, *popt)

plt.plot(xdata, ydata, 'o', label='实际数据')
plt.plot(xdata, fitted, label='拟合曲线')
plt.xlabel('时间')
plt.ylabel('感染比例')
plt.legend()
plt.show()

基于SIR模型的预测方法 #

当模型拟合效果达标后，可按以下步骤进行预测：

1. 确定预测时间范围：比如要预测未来10天，生成新的时间序列x_pred = np.arange(22, 32)（或者包含原有时间点的更长序列）。
2. 使用拟合参数计算预测值：用拟合得到的beta和gamma（即popt中的值）代入fit_odeint函数，计算预测的感染比例。
3. 更新初始条件（可选）：如果有最新的实际数据，可以用最新的S、I、R比例替换初始值，提升预测准确性。

示例预测代码：

# 预测未来10天
x_pred = np.arange(0, 32)  # 包含原有22天+新增10天
predicted = fit_odeint(x_pred, *popt)

plt.plot(xdata, ydata, 'o', label='实际数据')
plt.plot(x_pred, predicted, label='预测曲线', linestyle='--')
plt.xlabel('时间')
plt.ylabel('感染比例')
plt.legend()
plt.show()

注意：如果你的数据本身不符合经典SIR模型的假设（如存在重复感染、干预措施），即使修正代码，拟合效果也可能不理想，此时需要更换更复杂的模型（如SIRS、SEIR或带时变参数的SIR）。

内容的提问来源于stack exchange，提问作者Leshui