简介这份资源围绕传染病动力学建模中的随机微分方程展开面向流行病建模、公共卫生干预评估方向的研究人员及高校研究生帮助读者掌握环境噪声下疫情动态的数学刻画与数据驱动参数估计方法。压缩包共1个docx文件约8.09MB以文档形式完整呈现模型推导、算法流程与数值仿真细节。内容从随机微积分出发用标准布朗运动描述不可预见变化构建能更精确反映真实疫情特征的传播模型数据取自中国官方统计资料覆盖2021年初至2022年底共21个月的确诊、治愈与疫苗接种数据。方法上采用Metropolis-Hastings算法完成贝叶斯参数推断与最大似然估计并通过大样本抽样验证解的稳定性数值仿真环节给出Milstein格式等关键实现细节。文档还包含离散化模型公式、参数取值、初始数据表与方差最小化估计步骤便于读者复现建模与估计全过程。目前已有151人学习。1. 传染病动力学建模里的随机微分方程从确定性 SIR 到参数估计的完整落地路径确定性 SIR 模型跑出来的曲线太干净了干净到不真实。真实疫情数据里新增病例从来不是平滑下降的它有平台期、有反弹、有突然的尖峰。如果你只用常微分方程ODE去拟合残差图会告诉你一个残酷的事实模型系统性地低估了波动。随机微分方程SDE就是来解决这个问题的——它把随机性显式地写进微分方程让模型本身就能产生波动而不是把波动当成噪声扔掉。这篇笔记面向的是已经会写 SIR/SEIR 基础代码、想进一步做参数估计和不确定性量化的从业者。我会从 SDE 的构造方式讲起落到用真实疫情数据做参数估计的完整流程包括代码、参数设置和踩过的坑。核心问题只有一个怎么让传染病模型既保留动力学结构又能匹配真实数据的统计特征。2. 为什么要在传染病模型里引入随机项三种噪声来源与 SDE 构造方式2.1 确定性 SIR 的局限与噪声的三种来源确定性 SIR 的标准形式是 dS/dt -βSI/NdI/dt βSI/N - γIdR/dt γI。这套方程假设人群充分混合、参数恒定、没有外部干扰。但真实传播过程至少有三类随机性第一类是人口学随机性demographic stochasticity来源于个体接触和康复事件的离散性。当感染人数很少时每次接触是否导致传播本质上是一个伯努利试验感染人数少的时候这种随机性尤其显著可能直接决定疫情是否消亡。第二类是环境随机性environmental stochasticity来源于传播率 β 本身随时间波动。气候、人群行为变化、防控措施的松紧调整都会让 β 不是一个常数。这类随机性影响的是参数本身而不是事件计数。第三类是观测噪声来源于检测能力、报告延迟、无症状感染者的漏报。这类噪声不改变真实传播过程但会严重扭曲你看到的数据。SDE 主要处理的是前两类。常见做法是把环境随机性建模为对 β 的扰动β(t) β₀ σ·ξ(t)其中 ξ(t) 是白噪声。代入 SIR 后得到 Ito 型 SDEdS -β₀SI/N dt - σSI/N dW dI (β₀SI/N - γI) dt σSI/N dW dR γI dt这里 dW 是维纳过程增量。注意 σ 项同时出现在 S 和 I 的方程里因为传播事件同时消耗易感者和增加感染者两个方程必须共享同一个噪声源否则总人口 N 不守恒。2.2 Ito 与 Stratonovich选哪个不是玄学SDE 的数值求解和参数估计都依赖于你对随机积分的解释。Ito 积分和 Stratonovich 积分的区别在于被积函数在区间端点取值的方式。对于传染病模型我一般用 Ito 形式原因是Ito 形式下噪声项的期望为零E[dW] 0这使得矩方程更容易推导参数估计时的似然函数也更直接。如果你从连续时间的物理过程出发建模Stratonovich 可能更自然但如果你是从离散事件近似出发传染病传播本质上是离散事件Ito 是更合理的选择。两者之间可以通过一个修正项互相转换但对于 SIR 这类模型转换项会引入额外的漂移修正反而让参数解释变得复杂。所以我的建议是统一用 Ito不要混用。2.3 用 Python 构造 SIR-SDE 并做一次数值模拟下面是一个最小可运行的 SIR-SDE 模拟代码使用 Euler-Maruyama 方法import numpy as np import matplotlib.pyplot as plt def sir_sde_simulate(S0, I0, R0, beta, gamma, sigma, T, dt, seed42): Euler-Maruyama 方法模拟 SIR-SDE S0, I0, R0: 初始人数 beta: 基础传播率 gamma: 康复率 sigma: 噪声强度 T: 总模拟时间天 dt: 时间步长 rng np.random.default_rng(seed) N S0 I0 R0 n_steps int(T / dt) S np.zeros(n_steps 1) I np.zeros(n_steps 1) R np.zeros(n_steps 1) S[0], I[0], R[0] S0, I0, R0 for t in range(n_steps): dW rng.normal(0, np.sqrt(dt)) # 维纳过程增量 # 传播项 infection beta * S[t] * I[t] / N recovery gamma * I[t] # 噪声项同时作用于 S 和 I noise sigma * S[t] * I[t] / N * dW S[t1] S[t] - infection * dt - noise I[t1] I[t] (infection - recovery) * dt noise R[t1] R[t] recovery * dt # 边界保护人数不能为负 S[t1] max(S[t1], 0) I[t1] max(I[t1], 0) R[t1] max(R[t1], 0) return S, I, R # 参数设置 S0, I0, R0 9900, 100, 0 beta, gamma 0.3, 0.1 sigma 0.05 T, dt 160, 0.1 S, I, R sir_sde_simulate(S0, I0, R0, beta, gamma, sigma, T, dt)这段代码的关键点有三个。第一dW的标准差是sqrt(dt)不是dt这是维纳过程的定义决定的写错的话噪声强度会随步长变化结果不可复现。第二噪声项noise同时从 S 中减去并加到 I 中保证 SIR 守恒。第三边界保护max(..., 0)是必须的Euler-Maruyama 在步长较大时可能产生负值虽然理论上 SDE 的解不会为负但数值离散会。参数方面sigma控制噪声强度一般取 0.01 到 0.1 之间。dt建议不超过 0.1否则离散误差会显著影响结果。如果你要做参数估计dt需要和实际数据的时间粒度匹配比如数据是每日报告dt就设为 1。3. 基于疫情数据的参数估计从最小二乘到粒子滤波的完整实现3.1 参数估计的目标你到底在估什么在 SDE 框架下做参数估计需要区分三类参数第一类是动力学参数比如 β、γ、σ。这些参数决定了模型的动力学行为和噪声强度。第二类是初始条件S₀、I₀、R₀。实际数据中初始感染人数往往是未知的需要一起估计。第三类是观测参数比如报告率 ρ真实感染中被报告的比例、观测噪声方差。这些参数连接模型状态和实际观测。常见做法是固定一部分参数比如 γ 可以从临床数据估计潜伏期可以从文献获取只估计 β、σ 和初始条件。这样做的原因是同时估计所有参数会导致可识别性问题——不同的参数组合可能产生几乎相同的输出曲线。3.2 用最小二乘做初步估计快速但不够最小二乘是最直观的方法找到一组参数使得模型输出和观测数据的残差平方和最小。对于 SDE由于每次模拟的随机路径不同需要做多次模拟取平均或者用矩匹配的方法。from scipy.optimize import minimize import numpy as np def loss_function(params, data, n_sims50): params: [beta, gamma, sigma, I0] data: 观测的每日新增病例数组 n_sims: 每次评估时的模拟次数 beta, gamma, sigma, I0 params S0 100000 - I0 R0 0 T len(data) dt 1.0 sims [] for i in range(n_sims): _, I_sim, _ sir_sde_simulate(S0, I0, R0, beta, gamma, sigma, T, dt, seedi) # 计算每日新增I[t1] - I[t] gamma*I[t] 近似 daily_new np.diff(I_sim) gamma * I_sim[:-1] sims.append(daily_new[:len(data)]) sim_mean np.mean(sims, axis0) residual data - sim_mean return np.sum(residual**2) # 假设 data 是观测的每日新增病例 # result minimize(loss_function, x0[0.3, 0.1, 0.05, 100], # args(data,), methodNelder-Mead, # options{maxiter: 500})这段代码的逻辑是对每组候选参数模拟多条路径取平均后计算与观测数据的残差平方和。n_sims取 50 是一个折中太少会导致损失函数本身噪声太大优化器无法收敛太多则计算成本高。methodNelder-Mead适合这种目标函数没有解析梯度的情况。但最小二乘有两个问题。第一它假设观测噪声是高斯分布且方差恒定而实际疫情数据的噪声往往是非高斯的、异方差的。第二它只给出了点估计没有不确定性量化。你无法知道估计出来的 β 有多可靠。3.3 粒子滤波做似然估计给参数一个置信区间要解决不确定性问题需要从贝叶斯角度出发计算参数的后验分布。粒子滤波Particle Filter是处理非线性、非高斯状态空间模型的标准工具。核心思路是把 SDE 离散化后看作状态空间模型用一组粒子样本近似状态的后验分布然后通过观测数据不断更新权重最终得到参数的似然函数值。def particle_filter_loglik(params, data, n_particles1000, seed0): 用粒子滤波计算给定参数下的对数似然 params: [beta, gamma, sigma, I0] data: 观测的每日新增病例 rng np.random.default_rng(seed) beta, gamma, sigma, I0 params N 100000 S0 N - I0 R0 0 dt 1.0 # 初始化粒子 particles_S np.full(n_particles, S0, dtypefloat) particles_I np.full(n_particles, I0, dtypefloat) particles_R np.full(n_particles, R0, dtypefloat) log_weights np.zeros(n_particles) loglik 0.0 obs_noise_std 10.0 # 观测噪声标准差需要根据数据调整 for t in range(len(data)): # 传播步骤 dW rng.normal(0, np.sqrt(dt), n_particles) infection beta * particles_S * particles_I / N recovery gamma * particles_I noise sigma * particles_S * particles_I / N * dW particles_S np.maximum(particles_S - infection*dt - noise, 0) particles_I np.maximum(particles_I (infection - recovery)*dt noise, 0) particles_R np.maximum(particles_R recovery*dt, 0) # 计算预测的每日新增 pred_new infection * dt noise # 近似 # 更新权重观测模型为高斯 log_inc -0.5 * ((data[t] - pred_new) / obs_noise_std)**2 log_weights log_inc # 归一化并计算似然贡献 max_lw np.max(log_weights) weights np.exp(log_weights - max_lw) weights / np.sum(weights) loglik max_lw np.log(np.sum(np.exp(log_weights - max_lw))) - np.log(n_particles) # 重采样当有效粒子数过低时 ess 1.0 / np.sum(weights**2) if ess n_particles / 2: indices rng.choice(n_particles, sizen_particles, pweights) particles_S particles_S[indices] particles_I particles_I[indices] particles_R particles_R[indices] log_weights np.zeros(n_particles) return loglik这段代码的关键参数和逻辑说明n_particles取 1000 是一个经验值。粒子数太少似然估计的方差大粒子数太多计算成本高。对于参数估计通常需要 500 到 2000 个粒子。obs_noise_std是观测噪声的标准差需要根据数据的报告质量来设定。如果数据波动大这个值应该设大一些如果数据平滑设小一些。这个参数本身也可以作为待估参数但会增加计算复杂度。重采样步骤是粒子滤波的核心。当有效粒子数ESS低于阈值时根据权重重新采样避免粒子退化。阈值一般取n_particles / 2。得到对数似然后可以用 MCMC 采样参数的后验分布def mcmc_sampling(data, n_iter5000, burn_in1000): 简单的随机游走 MCMC rng np.random.default_rng(123) # 初始参数 [beta, gamma, sigma, I0] current np.array([0.3, 0.1, 0.05, 100.0]) current_loglik particle_filter_loglik(current, data) samples [] accepted 0 proposal_scale np.array([0.01, 0.005, 0.005, 5.0]) for i in range(n_iter): proposal current rng.normal(0, proposal_scale) # 参数约束 if proposal[0] 0 or proposal[1] 0 or proposal[2] 0 or proposal[3] 0: samples.append(current.copy()) continue proposal_loglik particle_filter_loglik(proposal, data) log_alpha proposal_loglik - current_loglik if np.log(rng.uniform()) log_alpha: current proposal current_loglik proposal_loglik accepted 1 samples.append(current.copy()) samples np.array(samples[burn_in:]) print(f接受率: {accepted / n_iter:.3f}) return samplesproposal_scale控制每步的跳跃幅度。太大接受率低太小混合慢。目标接受率在 0.2 到 0.4 之间。burn_in是预热期丢弃前面的样本让链收敛到后验分布。4. 避坑指南SDE 参数估计中最容易翻车的五个地方4.1 噪声项写错导致总人口不守恒现象模拟结果中 SIR 随时间漂移有时增加有时减少。原因噪声项没有同时作用于 S 和 I或者符号写反了。在 SIR-SDE 中传播事件同时减少 S 和增加 I噪声项必须共享同一个 dW且符号相反。解决检查代码中噪声项的符号。S 的方程中噪声项为负I 的方程中噪声项为正且两者使用同一个 dW 样本。R 的方程中不应有噪声项康复过程不引入额外的环境随机性除非你显式建模康复率的随机性。4.2 步长 dt 太大导致数值不稳定现象模拟结果出现负值或者 I(t) 突然爆炸。原因Euler-Maruyama 方法的收敛阶是 0.5强收敛步长太大时离散误差会累积。特别是当 sigma 较大时噪声项的贡献可能超过漂移项导致状态跳出合理范围。解决减小 dt一般建议 dt ≤ 0.1。如果计算成本不允许可以用 Milstein 方法收敛阶提高到 1.0但需要计算噪声项的导数。另一个实用技巧是自适应步长当状态接近边界时自动减小 dt。4.3 粒子滤波中粒子退化导致似然估计有偏现象MCMC 链不收敛或者后验分布异常窄。原因粒子滤波中权重集中在少数粒子上有效粒子数迅速下降。如果重采样阈值设得太低或者观测噪声设得太小都会加剧退化。解决监控 ESS确保重采样及时触发。观测噪声不要设得太小一般不低于数据标准差的 10%。如果仍然退化增加粒子数或改用辅助粒子滤波。4.4 参数可识别性问题导致估计结果无意义现象MCMC 链在不同参数值之间跳跃后验分布出现多峰或非常宽。原因β 和 σ 在某些数据模式下可能高度相关不同的 (β, σ) 组合产生相似的输出。初始条件 I0 和报告率 ρ 也容易混淆。解决固定一部分参数。γ 可以从临床数据估计报告率可以从血清学调查获取。如果必须估计所有参数需要加入强先验信息或者使用多组数据比如不同地区的疫情数据联合估计。4.5 观测模型与实际数据生成过程不匹配现象拟合效果看起来不错但预测区间完全覆盖不了未来的实际数据。原因观测模型假设高斯噪声但实际数据可能是过度离散的方差大于均值或者存在报告延迟和周期性周末效应。解决改用负二项观测模型处理过度离散加入报告延迟的卷积模型或者在观测模型中显式建模周末效应。这些改动会增加计算复杂度但能显著改善预测的可靠性。5. 进阶技巧用矩匹配加速参数估计并验证模型5.1 矩匹配绕开粒子滤波的快速近似粒子滤波虽然通用但计算成本高。如果你只需要点估计可以用矩匹配方法推导 SDE 的均值和方差的近似演化方程然后匹配观测数据的一阶和二阶矩。对于 SIR-SDE均值的演化近似为确定性 SIR 方程方差的演化可以通过线性化得到。具体做法是令 x [S, I, R]ᵀSDE 可以写成 dx f(x)dt g(x)dW。在均值附近做一阶泰勒展开得到均值和协方差的 ODEdμ/dt f(μ) dP/dt J(μ)P P J(μ)ᵀ g(μ)g(μ)ᵀ其中 J 是 f 的雅可比矩阵。这套方程可以用标准 ODE 求解器快速积分计算成本远低于粒子滤波。from scipy.integrate import odeint def moment_equations(y, t, beta, gamma, sigma, N): y [S_mean, I_mean, R_mean, P11, P12, P13, P22, P23, P33] 协方差矩阵对称只存上三角 S, I, R y[0], y[1], y[2] P11, P12, P13, P22, P23, P33 y[3:] # 均值方程 dS -beta * S * I / N dI beta * S * I / N - gamma * I dR gamma * I # 雅可比矩阵 J np.array([ [-beta*I/N, -beta*S/N, 0], [beta*I/N, beta*S/N - gamma, 0], [0, gamma, 0] ]) # 噪声矩阵 g(x)g(x)^T g np.array([-sigma*S*I/N, sigma*S*I/N, 0]) G np.outer(g, g) P np.array([[P11, P12, P13], [P12, P22, P23], [P13, P23, P33]]) dP J P P J.T G return [dS, dI, dR, dP[0,0], dP[0,1], dP[0,2], dP[1,1], dP[1,2], dP[2,2]]矩匹配的优点是快缺点是只用了前两阶矩如果数据分布严重偏离高斯估计会有偏。我一般用它做初步探索确定参数的大致范围再用粒子滤波做精细估计。5.2 后验预测检验你的模型真的靠谱吗参数估计完成后必须做后验预测检验。具体做法是从后验分布中抽取参数样本对每组参数模拟多条路径然后检查模拟数据是否能够复现观测数据的关键统计特征。我通常检查三个指标每日新增的均值和方差、疫情峰值时间和峰值大小、以及自相关函数。如果模拟数据的方差系统性低于观测数据说明噪声强度 σ 估计偏低如果峰值时间偏差大说明 β 或 γ 有偏。def posterior_predictive_check(samples, data, n_draws200): 从后验样本中抽取参数模拟并计算统计量 rng np.random.default_rng(456) indices rng.choice(len(samples), sizen_draws, replaceFalse) sim_stats [] for idx in indices: beta, gamma, sigma, I0 samples[idx] S0 100000 - I0 _, I_sim, _ sir_sde_simulate(S0, I0, 0, beta, gamma, sigma, len(data), 1.0, seedidx) daily_new np.diff(I_sim) gamma * I_sim[:-1] sim_stats.append([ np.mean(daily_new[:len(data)]), np.var(daily_new[:len(data)]), np.argmax(daily_new[:len(data)]) ]) sim_stats np.array(sim_stats) obs_stats [np.mean(data), np.var(data), np.argmax(data)] print(观测统计量:, obs_stats) print(模拟均值:, np.mean(sim_stats, axis0)) print(模拟95%区间:, np.percentile(sim_stats, [2.5, 97.5], axis0))如果观测统计量落在模拟的 95% 区间内说明模型没有系统性偏差。如果落在区间外需要检查模型设定是否有问题。5.3 一个我反复使用的习惯每次做参数估计之前我会先用模拟数据跑一遍完整流程给定已知参数生成模拟观测数据然后跑参数估计看能否恢复真实参数。这一步能帮你发现代码 bug、可识别性问题和先验设置问题。如果连模拟数据都恢复不了真实数据就更不用想了。这个习惯帮我省了很多时间。有一次我发现 β 的估计值总是偏高后来用模拟数据测试才发现是观测模型中报告率的先验设得太窄导致 β 被迫补偿。修正先验后估计就正常了。希望帮到你。本文还有配套的精品资源点击获取
