从Logistic曲线到疫情预测：用Python和SciPy复现SI传染病模型（附代码）

从Logistic曲线到疫情预测用Python和SciPy复现SI传染病模型附代码最近在整理疫情数据时我发现一个有趣的现象很多地区的感染人数增长曲线都呈现出典型的S型特征。这让我想起了经典的SI传染病模型它用简单的微分方程就能描述这种增长模式。今天我们就抛开复杂的数学推导直接用Python代码来实现这个模型看看如何用它来预测疫情发展趋势。1. 理解SI模型的核心思想SI模型是传染病建模中最基础的模型之一它把人群分为两类易感者(Susceptible)可能被感染的健康人群感染者(Infective)已经患病并具有传染性的人群模型的核心假设包括总人口数量恒定不考虑出生、死亡和迁移感染者无法被治愈或获得免疫每个感染者每天接触λ个其他人其中健康者会被感染注意虽然这些假设在现实中并不完全成立但简化后的模型仍能提供有价值的预测参考。根据这些假设我们可以建立微分方程来描述感染比例i(t)的变化di/dt λ * i(t) * (1 - i(t))这个方程的解就是著名的Logistic函数i(t) 1 / (1 (1/i0 - 1) * e^(-λt))其中i0是初始感染比例。这个函数呈现出典型的S型增长特征与很多实际疫情数据高度吻合。2. 搭建Python计算环境在开始编码前我们需要准备以下工具# 必需库安装如果尚未安装 # pip install numpy scipy matplotlib主要使用的库及其作用NumPy提供高效的数值计算支持SciPy用于求解微分方程Matplotlib数据可视化建议使用Jupyter Notebook进行交互式开发可以实时查看计算结果和图表。3. 数值求解SI模型微分方程虽然SI模型有解析解但为了后续更复杂模型的扩展我们先学习用数值方法求解。SciPy的odeint函数非常适合这类常微分方程问题。import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def si_model(i, t, lambda_): 定义SI模型的微分方程 return lambda_ * i * (1 - i) # 参数设置 lambda_ 0.3 # 日接触率 i0 0.01 # 初始感染比例 t np.linspace(0, 50, 100) # 时间范围0到50天 # 求解微分方程 solution odeint(si_model, i0, t, args(lambda_,)) i solution[:, 0]这段代码的核心是si_model函数它定义了微分方程的右侧。odeint会自动处理时间步进和数值积分返回每个时间点的感染比例。4. 可视化结果与分析让我们将数值解与理论解进行对比并分析关键特征# 计算理论解Logistic函数 theory_i 1 / (1 (1/i0 - 1) * np.exp(-lambda_ * t)) # 绘制结果 plt.figure(figsize(10, 6)) plt.plot(t, i, b-, linewidth2, label数值解) plt.plot(t, theory_i, r--, linewidth2, label理论解) plt.xlabel(时间(天)) plt.ylabel(感染比例) plt.title(SI模型感染曲线 (λ0.3)) plt.legend() plt.grid(True) # 标记关键点 peak_t np.log(1/i0 - 1)/lambda_ plt.axvline(peak_t, colorgray, linestyle:) plt.text(peak_t1, 0.2, f高峰期 t{peak_t:.1f}天, rotation90) plt.show()图表会显示两条几乎重合的曲线验证了我们数值解的正确性。关键观察点包括感染高峰期当感染比例达到50%时新增感染速度最快饱和期随着易感人群减少新增感染逐渐放缓最终状态所有人都被感染i→15. 参数λ对疫情发展的影响日接触率λ是模型中最关键的参数它反映了防控措施的严格程度。让我们比较不同λ值下的曲线变化lambda_values [0.1, 0.3, 0.5, 0.8] plt.figure(figsize(10, 6)) for lambda_ in lambda_values: solution odeint(si_model, i0, t, args(lambda_,)) plt.plot(t, solution[:, 0], labelfλ{lambda_}) plt.xlabel(时间(天)) plt.ylabel(感染比例) plt.title(不同日接触率下的感染曲线) plt.legend() plt.grid(True) plt.show()从图中可以直观看出λ值高峰期到来时间曲线陡峭程度现实对应措施0.1较晚平缓严格封锁0.3中等适中部分限制0.5较早较陡轻度防控0.8很早非常陡无防控措施这个分析告诉我们即使简单的SI模型也能为防控策略提供量化参考。通过社交距离等措施降低λ值可以有效延缓疫情高峰的到来。6. 模型局限性讨论虽然SI模型简单实用但在实际应用中需要注意以下限制无法考虑康复和免疫更适合描述无法治愈的传染病忽略潜伏期从暴露到具有传染性的过程未被建模同质混合假设现实中接触网络往往是非均匀的忽略人口动态变化出生、死亡和迁移未被考虑对于更复杂的情况可以考虑以下扩展方向SIR模型加入康复者群体SEIR模型增加潜伏期人群网络模型考虑接触网络结构7. 实战应用预测疫情发展让我们用一个实际案例演示如何使用SI模型进行预测。假设某地区初始感染率为1%通过早期数据估计λ≈0.25# 预测未来60天发展 lambda_ 0.25 i0 0.01 t_pred np.linspace(0, 60, 100) solution odeint(si_model, i0, t_pred, args(lambda_,)) i_pred solution[:, 0] # 计算关键指标 peak_t np.log(1/i0 - 1)/lambda_ peak_i 0.5 max_rate lambda_ * peak_i * (1 - peak_i) print(f预测结果 - 疫情高峰将在{peak_t:.1f}天后到来 - 最大日新增感染比例{max_rate:.2%})这个简单预测可以给决策者提供以下参考预计疫情高峰到来的时间点医疗资源需求的峰值时间不同防控措施可能带来的影响在实际项目中还需要结合更多数据定期更新λ的估计值以提高预测准确性。