用Python和MATLAB实战传染病模型：从SI到SEIR的疫情预测全解析
当新冠疫情席卷全球时，传染病模型突然从学术论文走进了公众视野。作为数据科学家或开发者，我们不仅要理解这些模型背后的数学原理，更需要掌握如何用代码将它们转化为可交互的模拟工具。本文将带你用Python和MATLAB两种主流工具，完整实现五种经典传染病模型，并通过可视化分析不同干预措施的效果。
1. 环境准备与工具选择
在开始建模前，我们需要配置合适的开发环境。Python和MATLAB各有优势：Python开源免费且生态丰富，MATLAB则在数学建模和仿真方面有着深厚的积累。
Python环境配置：
# 创建虚拟环境 python -m venv epi-model source epi-model/bin/activate # Linux/Mac epi-model\Scripts\activate # Windows # 安装必要库 pip install numpy scipy matplotlib pandas ipythonMATLAB准备：
- 确保已安装Symbolic Math Toolbox和Statistics and Machine Learning Toolbox
- 推荐使用Live Script格式，便于交互式开发和文档记录
提示：对于教育用户，MATLAB Online提供免费的云端计算环境，适合没有本地授权的学习者
两种工具的核心差异对比如下：
| 特性 | Python | MATLAB |
|---|---|---|
| 微分方程求解 | scipy.integrate.odeint | ode45/ode23系列 |
| 符号计算 | sympy库（速度较慢） | 内置Symbolic Math Toolbox |
| 可视化 | matplotlib/seaborn | 内置绘图函数 |
| 性能 | 依赖实现方式 | 矩阵运算优化好 |
| 学习曲线 | 需组合多个库 | 一体化体验 |
2. 基础模型实现：从SI到SIR
2.1 SI模型：最简单的传染场景
SI模型假设人群只分为易感者(S)和感染者(I)，没有恢复和免疫机制。这虽然简单，却是理解传染病动力学的基础。
Python实现：
import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def si_model(y, t, beta, N): S, I = y dSdt = -beta * S * I / N dIdt = beta * S * I / N return [dSdt, dIdt] # 参数设置 N = 1000 # 总人口 I0 = 1 # 初始感染者 beta = 0.3 # 传染率 t = np.linspace(0, 100, 100) # 100天模拟 # 求解微分方程 solution = odeint(si_model, [N-I0, I0], t, args=(beta, N)) S, I = solution.T # 可视化 plt.figure(figsize=(10,6)) plt.plot(t, S, label='Susceptible') plt.plot(t, I, label='Infected') plt.xlabel('Days') plt.ylabel('Population') plt.legend() plt.title('SI Model Simulation') plt.grid(True)MATLAB实现：
function si_model_demo % 参数设置 N = 1000; % 总人口 I0 = 1; % 初始感染者 beta = 0.3; % 传染率 tspan = 0:1:100; % 时间范围 % 定义微分方程 si_ode = @(t,y) [-beta*y(1)*y(2)/N; beta*y(1)*y(2)/N]; % 求解 [t,y] = ode45(si_ode, tspan, [N-I0; I0]); % 绘图 figure plot(t, y(:,1), 'b-', t, y(:,2), 'r-') xlabel('Days') ylabel('Population') legend('Susceptible','Infected') title('SI Model Simulation') grid on end2.2 SIR模型：引入康复机制
SIR模型在SI基础上增加了康复者(R)群体，更接近真实传染病的发展规律。我们需要新增康复率γ参数。
关键参数解释：
- 基本传染数R0 = β/γ，表示一个感染者平均传染多少人
- 当R0 > 1时，疫情会扩散；R0 < 1时，疫情会逐渐消失
Python代码升级：
def sir_model(y, t, beta, gamma, N): S, I, R = y dSdt = -beta * S * I / N dIdt = beta * S * I / N - gamma * I dRdt = gamma * I return [dSdt, dIdt, dRdt] # 新增gamma参数 gamma = 0.05 solution = odeint(sir_model, [N-I0, I0, 0], t, args=(beta, gamma, N)) S, I, R = solution.T3. 进阶模型：应对复杂传染场景
3.1 SEIR模型：考虑潜伏期
新冠疫情让我们认识到潜伏期(E)的重要性。SEIR模型在SIR基础上增加了暴露人群，更准确地模拟了类似COVID-19的传染病。
模型参数扩展：
- σ：潜伏期转化为感染者的速率（σ = 1/平均潜伏期）
- 新增潜伏期人群E的动态变化
MATLAB实现要点：
function seir_model % 参数设置 N = 1e6; beta = 0.5; sigma = 1/5.2; gamma = 1/12; tspan = 0:1:180; seir_ode = @(t,y) [ -beta*y(1)*y(3)/N; beta*y(1)*y(3)/N - sigma*y(2); sigma*y(2) - gamma*y(3); gamma*y(3) ]; [t,y] = ode45(seir_ode, tspan, [N-10 0 10 0]); % 可视化代码... end3.2 带干预措施的SIRS模型
现实中的防疫措施会改变模型参数。我们可以通过时间依赖的参数模拟封控、社交隔离等干预：
def time_varying_beta(t): """模拟不同阶段的防控措施""" if t < 30: return 0.5 # 正常传播 elif t < 60: return 0.2 # 严格封控 else: return 0.3 # 部分解封 def sirs_model(y, t, gamma, delta, N): S, I, R = y beta = time_varying_beta(t) dSdt = -beta * S * I / N + delta * R dIdt = beta * S * I / N - gamma * I dRdt = gamma * I - delta * R return [dSdt, dIdt, dRdt]4. 模型验证与参数估计
构建模型只是第一步，如何确定参数和验证模型准确性同样重要。我们可以使用真实数据进行参数拟合。
Python参数估计示例：
from scipy.optimize import minimize def fit_sir(params, infected_data, N): beta, gamma = params def sir_ode(y, t): ... # 同上 solution = odeint(sir_ode, [N-1, 1, 0], range(len(infected_data)), args=(beta, gamma, N)) predicted = solution[:,1] return np.sum((predicted - infected_data)**2) # MSE # 使用某地区真实感染数据 data = np.array([...]) initial_guess = [0.3, 0.1] result = minimize(fit_sir, initial_guess, args=(data, 1e6)) beta_opt, gamma_opt = result.x关键验证指标：
- 均方误差(MSE)：衡量模型预测与真实数据的差距
- R0估计：评估传染性强弱
- 峰值时间预测：判断疫情拐点
5. 高级应用与可视化技巧
5.1 交互式疫情模拟
使用Python的ipywidgets创建可调节参数的交互式模拟：
from ipywidgets import interact def plot_sir(beta=0.3, gamma=0.1, N=1000, days=100): t = np.linspace(0, days, days) solution = odeint(sir_model, [N-1, 1, 0], t, args=(beta, gamma, N)) # 绘图代码... interact(plot_sir, beta=(0.01, 1.0, 0.01), gamma=(0.01, 0.5, 0.01), N=(100, 10000, 100), days=(30, 365, 10))5.2 三维参数空间探索
使用MATLAB可视化R0对疫情发展的影响：
[R0_grid, gamma_grid] = meshgrid(0.5:0.1:3, 0.01:0.01:0.2); peak_infected = zeros(size(R0_grid)); for i = 1:numel(R0_grid) beta = R0_grid(i)*gamma_grid(i); [~,y] = ode45(@(t,y) sir_ode(t,y,beta,gamma_grid(i),1e6), 0:1:365, [1e6-100 100 0]); peak_infected(i) = max(y(:,2)); end figure surf(R0_grid, gamma_grid, peak_infected) xlabel('R0'); ylabel('\gamma'); zlabel('Peak Infections') title('Epidemic Peak vs R0 and Recovery Rate')5.3 地理信息集成
将模型结果与地图数据结合，展示疫情空间传播：
import geopandas as gpd def regional_simulation(regions): results = {} for region in regions: # 为每个地区设置不同参数 beta = region['density'] * 0.001 solution = odeint(sir_model, ...) results[region['name']] = solution # 在地图上可视化 world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) ax = world.plot(figsize=(15,10)) for region, data in results.items(): peak_day = np.argmax(data[:,1]) # 在地图上标记各地区的疫情峰值时间 # ...在实际项目中，我发现模型的准确性高度依赖于参数的选择。特别是在模拟COVID-19时，将社交隔离措施量化为β参数的动态变化是关键挑战。一个实用的技巧是使用移动平均处理实时感染数据，再将其输入到参数估计算法中，这样可以减少数据报告延迟带来的影响。
