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

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

最近在整理疫情数据时，我发现一个有趣的现象：很多地区的感染人数增长曲线都呈现出典型的S型特征。这让我想起了经典的SI传染病模型，它用简单的微分方程就能描述这种增长模式。今天，我们就抛开复杂的数学推导，直接用Python代码来实现这个模型，看看如何用它来预测疫情发展趋势。

1. 理解SI模型的核心思想

SI模型是传染病建模中最基础的模型之一，它把人群分为两类：

• 易感者(Susceptible)：可能被感染的健康人群
• 感染者(Infective)：已经患病并具有传染性的人群

模型的核心假设包括：

1. 总人口数量恒定（不考虑出生、死亡和迁移）
2. 感染者无法被治愈或获得免疫
3. 每个感染者每天接触λ个其他人，其中健康者会被感染

注意：虽然这些假设在现实中并不完全成立，但简化后的模型仍能提供有价值的预测参考。

根据这些假设，我们可以建立微分方程来描述感染比例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-', linewidth=2, label='数值解') plt.plot(t, theory_i, 'r--', linewidth=2, 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, color='gray', linestyle=':') plt.text(peak_t+1, 0.2, f'高峰期 t={peak_t:.1f}天', rotation=90) plt.show()

图表会显示两条几乎重合的曲线，验证了我们数值解的正确性。关键观察点包括：

• 感染高峰期：当感染比例达到50%时，新增感染速度最快
• 饱和期：随着易感人群减少，新增感染逐渐放缓
• 最终状态：所有人都被感染（i→1）

5. 参数λ对疫情发展的影响

日接触率λ是模型中最关键的参数，它反映了防控措施的严格程度。让我们比较不同λ值下的曲线变化：

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], label=f'λ={lambda_}') plt.xlabel('时间(天)') plt.ylabel('感染比例') plt.title('不同日接触率下的感染曲线') plt.legend() plt.grid(True) plt.show()

从图中可以直观看出：

这个分析告诉我们，即使简单的SI模型也能为防控策略提供量化参考。通过社交距离等措施降低λ值，可以有效延缓疫情高峰的到来。

6. 模型局限性讨论

虽然SI模型简单实用，但在实际应用中需要注意以下限制：

1. 无法考虑康复和免疫：更适合描述无法治愈的传染病
2. 忽略潜伏期：从暴露到具有传染性的过程未被建模
3. 同质混合假设：现实中接触网络往往是非均匀的
4. 忽略人口动态变化：出生、死亡和迁移未被考虑

对于更复杂的情况，可以考虑以下扩展方向：

• 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%}")

这个简单预测可以给决策者提供以下参考：

• 预计疫情高峰到来的时间点
• 医疗资源需求的峰值时间
• 不同防控措施可能带来的影响

在实际项目中，还需要结合更多数据定期更新λ的估计值，以提高预测准确性。