用Python实战传染病模型从SIR到SEIR的完整代码实现附可视化传染病模型是理解疾病传播动态的重要工具。在公共卫生领域这些模型帮助预测疫情发展趋势、评估防控措施效果。本文将带你用Python实现经典的SIR和SEIR模型通过代码演示从基础理论到实际应用的完整流程。1. 环境准备与基础概念在开始编码前我们需要配置Python环境并理解模型的核心参数。推荐使用Anaconda创建独立环境conda create -n epidemic python3.9 conda activate epidemic pip install numpy scipy matplotlib pandasSIR模型将人群分为三类S(Susceptible)易感人群I(Infectious)感染者R(Recovered)康复/免疫人群关键参数定义参数描述典型范围β传染率0.1-0.5γ康复率0.05-0.2R₀基本再生数β/γ提示R₀1表示疫情会扩散R₀1则疫情将逐渐消退2. SIR模型实现与可视化2.1 模型微分方程实现使用SciPy的odeint求解微分方程系统import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def sir_model(y, t, beta, gamma): S, I, R y dSdt -beta * S * I dIdt beta * S * I - gamma * I dRdt gamma * I return [dSdt, dIdt, dRdt]2.2 参数设置与求解配置初始条件和时间范围# 参数设置 N 10000 # 总人口 I0 1 # 初始感染者 beta 0.3 # 传染率 gamma 0.1 # 康复率 days 160 # 模拟天数 # 初始条件 S0 N - I0 R0 0 y0 [S0/N, I0/N, R0/N] # 标准化为比例 # 时间点 t np.linspace(0, days, days)求解并可视化结果# 求解微分方程 solution odeint(sir_model, y0, t, args(beta, gamma)) S, I, R solution.T * N # 还原为实际人数 # 可视化 plt.figure(figsize(10,6)) plt.plot(t, S, b, label易感者) plt.plot(t, I, r, label感染者) plt.plot(t, R, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(SIR模型传染病传播模拟) plt.legend() plt.grid() plt.show()典型输出曲线特征感染者数量呈现钟形曲线最终康复者比例取决于R₀值疫情持续时间与γ值负相关3. SEIR模型进阶实现3.1 引入潜伏期参数SEIR模型增加了暴露者(E)群体def seir_model(y, t, beta, gamma, sigma): S, E, I, R y dSdt -beta * S * I dEdt beta * S * I - sigma * E dIdt sigma * E - gamma * I dRdt gamma * I return [dSdt, dEdt, dIdt, dRdt]关键新增参数σ (sigma)潜伏期倒数1/σ为平均潜伏期典型取值σ0.2 (对应5天潜伏期)3.2 完整SEIR实现代码# 参数设置 sigma 0.2 # 潜伏率 E0 0 # 初始潜伏者 y0_seir [(N-I0)/N, E0/N, I0/N, R0/N] # 求解SEIR模型 solution_seir odeint(seir_model, y0_seir, t, args(beta, gamma, sigma)) S_seir, E_seir, I_seir, R_seir solution_seir.T * N # 可视化 plt.figure(figsize(10,6)) plt.plot(t, S_seir, b, label易感者) plt.plot(t, E_seir, m, label潜伏者) plt.plot(t, I_seir, r, label感染者) plt.plot(t, R_seir, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(SEIR模型传染病传播模拟) plt.legend() plt.grid() plt.show()SEIR与SIR模型对比差异疫情峰值出现时间延迟初期增长曲线更为平缓更适合模拟有显著潜伏期的疾病4. 参数敏感性分析与实际应用4.1 传染率β的影响beta_values [0.2, 0.3, 0.4] colors [r, g, b] plt.figure(figsize(10,6)) for beta, color in zip(beta_values, colors): solution odeint(sir_model, y0, t, args(beta, gamma)) I solution[:,1] * N plt.plot(t, I, color, labelfβ{beta}) plt.xlabel(天数) plt.ylabel(感染人数) plt.title(不同传染率下的感染曲线) plt.legend() plt.grid() plt.show()4.2 交互式参数探索使用ipywidgets创建交互界面from ipywidgets import interact def plot_sir_interactive(beta0.3, gamma0.1, sigma0.2): # SEIR模型求解 solution odeint(seir_model, y0_seir, t, args(beta, gamma, sigma)) S, E, I, R solution.T * N # 绘图 plt.figure(figsize(10,6)) plt.plot(t, S, b, label易感者) plt.plot(t, E, m, label潜伏者) plt.plot(t, I, r, label感染者) plt.plot(t, R, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(fSEIR模型 (β{beta}, γ{gamma}, σ{sigma})) plt.legend() plt.grid() plt.show() interact(plot_sir_interactive, beta(0.1, 0.5, 0.05), gamma(0.05, 0.3, 0.05), sigma(0.1, 0.5, 0.05))4.3 实际应用建议参数校准方法使用历史数据通过最小二乘法拟合考虑使用MCMC方法进行参数估计模型扩展方向添加年龄分层结构考虑空间异质性引入疫苗接种因素# 参数校准示例框架 from scipy.optimize import minimize def calibrate_params(params, real_data): beta, gamma params solution odeint(sir_model, y0, t, args(beta, gamma)) predicted_I solution[:,1] return np.sum((predicted_I - real_data)**2) # 假设real_data是实际观测的感染比例 initial_guess [0.3, 0.1] result minimize(calibrate_params, initial_guess, args(real_data,)) optimal_beta, optimal_gamma result.x5. 高级可视化与结果分析5.1 动态传播过程展示使用Matplotlib动画功能from matplotlib.animation import FuncAnimation fig, ax plt.subplots(figsize(10,6)) line_S, ax.plot([], [], b, label易感者) line_I, ax.plot([], [], r, label感染者) line_R, ax.plot([], [], g, label康复者) def init(): ax.set_xlim(0, days) ax.set_ylim(0, N) ax.set_xlabel(天数) ax.set_ylabel(人数) ax.set_title(SIR模型动态传播过程) ax.legend() ax.grid() return line_S, line_I, line_R def update(frame): line_S.set_data(t[:frame], S[:frame]) line_I.set_data(t[:frame], I[:frame]) line_R.set_data(t[:frame], R[:frame]) return line_S, line_I, line_R ani FuncAnimation(fig, update, frameslen(t), init_funcinit, blitTrue, interval50) plt.close()5.2 关键指标计算peak_day t[np.argmax(I)] peak_infections np.max(I) final_recovered R[-1] print(f疫情峰值出现在第{peak_day:.0f}天) print(f最高感染人数{peak_infections:.0f}人) print(f最终康复比例{final_recovered/N:.1%})5.3 三维参数空间探索from mpl_toolkits.mplot3d import Axes3D beta_range np.linspace(0.1, 0.5, 20) gamma_range np.linspace(0.05, 0.2, 20) peak_infections np.zeros((len(beta_range), len(gamma_range))) for i, beta in enumerate(beta_range): for j, gamma in enumerate(gamma_range): solution odeint(sir_model, y0, t, args(beta, gamma)) peak_infections[i,j] np.max(solution[:,1]) X, Y np.meshgrid(gamma_range, beta_range) fig plt.figure(figsize(10,7)) ax fig.add_subplot(111, projection3d) ax.plot_surface(X, Y, peak_infections, cmapviridis) ax.set_xlabel(康复率 γ) ax.set_ylabel(传染率 β) ax.set_zlabel(峰值感染比例) ax.set_title(参数空间对疫情峰值的影响) plt.show()
用Python实战传染病模型：从SIR到SEIR的完整代码实现（附可视化）
用Python实战传染病模型从SIR到SEIR的完整代码实现附可视化传染病模型是理解疾病传播动态的重要工具。在公共卫生领域这些模型帮助预测疫情发展趋势、评估防控措施效果。本文将带你用Python实现经典的SIR和SEIR模型通过代码演示从基础理论到实际应用的完整流程。1. 环境准备与基础概念在开始编码前我们需要配置Python环境并理解模型的核心参数。推荐使用Anaconda创建独立环境conda create -n epidemic python3.9 conda activate epidemic pip install numpy scipy matplotlib pandasSIR模型将人群分为三类S(Susceptible)易感人群I(Infectious)感染者R(Recovered)康复/免疫人群关键参数定义参数描述典型范围β传染率0.1-0.5γ康复率0.05-0.2R₀基本再生数β/γ提示R₀1表示疫情会扩散R₀1则疫情将逐渐消退2. SIR模型实现与可视化2.1 模型微分方程实现使用SciPy的odeint求解微分方程系统import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def sir_model(y, t, beta, gamma): S, I, R y dSdt -beta * S * I dIdt beta * S * I - gamma * I dRdt gamma * I return [dSdt, dIdt, dRdt]2.2 参数设置与求解配置初始条件和时间范围# 参数设置 N 10000 # 总人口 I0 1 # 初始感染者 beta 0.3 # 传染率 gamma 0.1 # 康复率 days 160 # 模拟天数 # 初始条件 S0 N - I0 R0 0 y0 [S0/N, I0/N, R0/N] # 标准化为比例 # 时间点 t np.linspace(0, days, days)求解并可视化结果# 求解微分方程 solution odeint(sir_model, y0, t, args(beta, gamma)) S, I, R solution.T * N # 还原为实际人数 # 可视化 plt.figure(figsize(10,6)) plt.plot(t, S, b, label易感者) plt.plot(t, I, r, label感染者) plt.plot(t, R, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(SIR模型传染病传播模拟) plt.legend() plt.grid() plt.show()典型输出曲线特征感染者数量呈现钟形曲线最终康复者比例取决于R₀值疫情持续时间与γ值负相关3. SEIR模型进阶实现3.1 引入潜伏期参数SEIR模型增加了暴露者(E)群体def seir_model(y, t, beta, gamma, sigma): S, E, I, R y dSdt -beta * S * I dEdt beta * S * I - sigma * E dIdt sigma * E - gamma * I dRdt gamma * I return [dSdt, dEdt, dIdt, dRdt]关键新增参数σ (sigma)潜伏期倒数1/σ为平均潜伏期典型取值σ0.2 (对应5天潜伏期)3.2 完整SEIR实现代码# 参数设置 sigma 0.2 # 潜伏率 E0 0 # 初始潜伏者 y0_seir [(N-I0)/N, E0/N, I0/N, R0/N] # 求解SEIR模型 solution_seir odeint(seir_model, y0_seir, t, args(beta, gamma, sigma)) S_seir, E_seir, I_seir, R_seir solution_seir.T * N # 可视化 plt.figure(figsize(10,6)) plt.plot(t, S_seir, b, label易感者) plt.plot(t, E_seir, m, label潜伏者) plt.plot(t, I_seir, r, label感染者) plt.plot(t, R_seir, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(SEIR模型传染病传播模拟) plt.legend() plt.grid() plt.show()SEIR与SIR模型对比差异疫情峰值出现时间延迟初期增长曲线更为平缓更适合模拟有显著潜伏期的疾病4. 参数敏感性分析与实际应用4.1 传染率β的影响beta_values [0.2, 0.3, 0.4] colors [r, g, b] plt.figure(figsize(10,6)) for beta, color in zip(beta_values, colors): solution odeint(sir_model, y0, t, args(beta, gamma)) I solution[:,1] * N plt.plot(t, I, color, labelfβ{beta}) plt.xlabel(天数) plt.ylabel(感染人数) plt.title(不同传染率下的感染曲线) plt.legend() plt.grid() plt.show()4.2 交互式参数探索使用ipywidgets创建交互界面from ipywidgets import interact def plot_sir_interactive(beta0.3, gamma0.1, sigma0.2): # SEIR模型求解 solution odeint(seir_model, y0_seir, t, args(beta, gamma, sigma)) S, E, I, R solution.T * N # 绘图 plt.figure(figsize(10,6)) plt.plot(t, S, b, label易感者) plt.plot(t, E, m, label潜伏者) plt.plot(t, I, r, label感染者) plt.plot(t, R, g, label康复者) plt.xlabel(天数) plt.ylabel(人数) plt.title(fSEIR模型 (β{beta}, γ{gamma}, σ{sigma})) plt.legend() plt.grid() plt.show() interact(plot_sir_interactive, beta(0.1, 0.5, 0.05), gamma(0.05, 0.3, 0.05), sigma(0.1, 0.5, 0.05))4.3 实际应用建议参数校准方法使用历史数据通过最小二乘法拟合考虑使用MCMC方法进行参数估计模型扩展方向添加年龄分层结构考虑空间异质性引入疫苗接种因素# 参数校准示例框架 from scipy.optimize import minimize def calibrate_params(params, real_data): beta, gamma params solution odeint(sir_model, y0, t, args(beta, gamma)) predicted_I solution[:,1] return np.sum((predicted_I - real_data)**2) # 假设real_data是实际观测的感染比例 initial_guess [0.3, 0.1] result minimize(calibrate_params, initial_guess, args(real_data,)) optimal_beta, optimal_gamma result.x5. 高级可视化与结果分析5.1 动态传播过程展示使用Matplotlib动画功能from matplotlib.animation import FuncAnimation fig, ax plt.subplots(figsize(10,6)) line_S, ax.plot([], [], b, label易感者) line_I, ax.plot([], [], r, label感染者) line_R, ax.plot([], [], g, label康复者) def init(): ax.set_xlim(0, days) ax.set_ylim(0, N) ax.set_xlabel(天数) ax.set_ylabel(人数) ax.set_title(SIR模型动态传播过程) ax.legend() ax.grid() return line_S, line_I, line_R def update(frame): line_S.set_data(t[:frame], S[:frame]) line_I.set_data(t[:frame], I[:frame]) line_R.set_data(t[:frame], R[:frame]) return line_S, line_I, line_R ani FuncAnimation(fig, update, frameslen(t), init_funcinit, blitTrue, interval50) plt.close()5.2 关键指标计算peak_day t[np.argmax(I)] peak_infections np.max(I) final_recovered R[-1] print(f疫情峰值出现在第{peak_day:.0f}天) print(f最高感染人数{peak_infections:.0f}人) print(f最终康复比例{final_recovered/N:.1%})5.3 三维参数空间探索from mpl_toolkits.mplot3d import Axes3D beta_range np.linspace(0.1, 0.5, 20) gamma_range np.linspace(0.05, 0.2, 20) peak_infections np.zeros((len(beta_range), len(gamma_range))) for i, beta in enumerate(beta_range): for j, gamma in enumerate(gamma_range): solution odeint(sir_model, y0, t, args(beta, gamma)) peak_infections[i,j] np.max(solution[:,1]) X, Y np.meshgrid(gamma_range, beta_range) fig plt.figure(figsize(10,7)) ax fig.add_subplot(111, projection3d) ax.plot_surface(X, Y, peak_infections, cmapviridis) ax.set_xlabel(康复率 γ) ax.set_ylabel(传染率 β) ax.set_zlabel(峰值感染比例) ax.set_title(参数空间对疫情峰值的影响) plt.show()
