简介围绕传染病动力学建模中的随机微分方程这份文档系统介绍了利用标准布朗运动刻画环境随机性的建模思路构建了易感-感染-恢复SIR类随机传染病模型并结合中国官方新冠疫情统计数据采用Metropolis-Hastings算法完成贝叶斯框架下的最大似然参数估计。内容涵盖模型推导、Milstein格式离散化、参数估计迭代步骤与数值仿真细节同时给出S、I、R变量的递推公式和方差最小化准则等关键环节适合流行病建模研究人员、公共卫生政策制定者及高校研究生学习参考。资源包共1个docx文档大小8.09MB内部包含完整的模型公式、参数表与计算过程便于读者复现实验。目前已有151人学习浏览文档在数据驱动建模与随机微积分应用方面具有较高参考价值可为预测传染病发展趋势和评估干预措施提供方法支撑。1. 传染病动力学建模里的随机微分方程这套资源到底解决什么问题新冠疫情数据摆在面前S、I、R 三类人群占比看起来有波动但标准 SIR 模型是确定性常微分方程给定初值和参数轨迹是唯一一条曲线根本解释不了数据里的来回震荡——比如同一季度占比相同下一季度却突然跳变。这套资源的核心思路是把环境扰动塞进模型用随机微分方程SDE描述传染过程再用 2021 年底到 2024 年三季度的真实疫情数据反推模型参数 α接触传染率、δ免疫丧失率、ω康复率。它适合两类人一类是做传染病建模的研究生想从确定性模型转向随机模型缺一份能跑通全流程的参考另一类是公共卫生或政策评估相关从业者需要根据真实数据做参数反演而不是只停留在理论推导。我拆完这份资源的结论是模型本身不复杂Milstein 离散化是数学上严密的桥梁真正的难度在参数估计——目标函数不平滑、局部极小多、参数之间相互补偿这才是值得花时间琢磨的地方。2. 为什么标准 SIR 不够用随机扰动进入动力学方程的完整推演2.1 确定性模型的局限数据震荡是噪声还是动力学行为标准 SIR 模型写作dS/dt b - αkS⟨k⟩⁻¹Σ p(j)I_j - dS δRdI/dt αkS⟨k⟩⁻¹Σ p(j)I_j - (ωd)IdR/dt ωI - (δd)R其中 bd0.1 是出生率和自然死亡率⟨k⟩1.3418 是网络平均度p(j) 是度分布相关项。这套确定性框架默认一个前提每天新增感染人数完全由 α、δ、ω 决定没有任何不可控因素。但真实疫情数据不是这样——表一里的 I_k 序列从 0.083 跳到 0.333 再回落这种波动幅度在确定性 SIR 里要改变参数才能产生问题是你不可能每个季度都改一次参数。随机扰动进入模型后状态变量变成随机过程每条模拟轨迹都是一次可能的疫情走向。这对应一个直观事实同等防控强度下感染人数也有随机涨落。模型里引入标准布朗运动 W(t)噪声强度由 σ₁k、σ₂k、σ₃k 控制分别作用于 S、I、R 三个变量。理论上这就把原来的常微分方程组改写成了 Itô 型随机微分方程组dS [b - αkS⟨k⟩⁻¹Σp(j)I_j - dS δR]dt σ₁k S dW₁dI [αkS⟨k⟩⁻¹Σp(j)I_j - (ωd)I]dt σ₂k I dW₂dR [ωI - (δd)R]dt σ₃k R dW₃扰动项的乘法结构噪声乘在 S、I、R 自身上不是随便选的。它保证状态变量为 0 时噪声项也为 0不会出现感染人数已经清零又被负噪声拉回正值的荒谬情形。2.2 接触核 p(j) 的归一化手算一遍才能避开的隐藏坑模型里 p(j)ck⁻ʳ满足 Σⱼ₌₁ⁿ p(j)1取 r3、n5。这一步看起来只是归一化实际上决定了多层网络接触的权重形状。我把手算过程拆开Σⱼ₌₁⁵ j⁻³ 1 1/8 1/27 1/64 1/125 ≈ 1 0.125 0.03704 0.015625 0.008 1.185665所以 c ≈ 0.8434。注意 Σⱼ₌₁ⁿ j p(j) 1 才是完整条件只做 Σp(j)1 会差一个均值约束。完整做法是先算 c 再做校验Σ j·cj⁻³ ≈ 0.8434 × (1 2/8 3/27 4/64 5/125) ≈ 0.8434 × (10.250.11110.06250.04) ≈ 0.8434 × 1.4636 ≈ 1.234再按归一化比例归一化。这个步骤在每个时间步都要用到建议单独写个函数返回 p(j) 数组别在迭代循环里重复算。2.3 状态变量的总量约束与参数取值边界模型要求 SIR1 恒成立参数取 bd0.1、σ₁k0.2、σ₂k0.1。这里有个容易被忽略的设计σ₁k 比 σ₂k 大一倍意味着易感人群的随机涨落更强。逻辑上说得通——易感人群基数大检测、隔离政策变化首先冲击的是 S 的波动。总量约束在确定性模型里是天然满足的因为三式相加恰好把转移项全部抵消。但引入随机项后三个独立的布朗运动各自乘以不同的 σ三项之和不再精确等于 1。处理方式是模拟 S 和 IR1-S-I 兜底或者在每个 Milstein 步结束时做归一化投影。下面第 3 章的实现里会给出具体做法。3. Milstein 离散化把随机微分方程变成能跑的迭代公式3.1 为什么不用欧拉-Maruyama收敛阶数不够导致的系统性偏差SDE 的数值解比 ODE 麻烦核心原因是 Itô 积分的二阶项不能直接丢弃。对 SDE 做泰勒-Itô 展开保留到一阶项得到欧拉-Maruyama 格式强收敛阶只有 0.5保留到二阶项得到 Milstein 格式强收敛阶提升到 1.0。在 σ0.2 的噪声水平下欧拉格式每步会引入与 Δt 同阶的偏差迭代几百步后累计误差会扭曲参数估计结果——你反推出来的 α 里混进了离散化误差而不是真实的传染率。Milstein 格式多出来的是一个与噪声导数相关的修正项数学形式上是 (σ²/2)·(Δt)(ν²-1) 那一坨作用是把布朗运动路径的曲率补回来。3.2 差分方程拆解与完整 Python 实现原始资源里给出的 Milstein 离散格式我把三个方程完整列出来S_{k,i1} S_{k,i} Δt[b - αkS⟨k⟩⁻¹Σⱼ₌₁ⁿ p(j)I_{j,i} - dS_{k,i} δR_{k,i}] σ₁k S_{k,i} Δt^{1/2} ν_{k,i} (σ₁k²/2)S_{k,i} Δt(ν_{k,i}² - 1)I_{k,i1} I_{k,i} Δt[αkS⟨k⟩⁻¹Σⱼ₌₁ⁿ p(j)I_{j,i} - (ωd)I_{k,i}] σ₂k I_{k,i} Δt^{1/2} ν_{k,i} (σ₂k²/2)I_{k,i} Δt(ν_{k,i}² - 1)R_{k,i1} R_{k,i} Δt[ωI_{k,i} - (δd)R_{k,i}] σ₃k R_{k,i} Δt^{1/2} ν_{k,i} (σ₃k²/2)R_{k,i} Δt(ν_{k,i}² - 1)其中 ν_{k,i} 是独立同分布的标准正态随机变量。注意公式里 σ³k 的取值资源里没有明确给出我用 σ₁k 和 σ₂k 的中间值 0.15 作为默认复现时你可以按数据波动程度调整。import numpy as np def compute_pj(n5, r3): 接触核 p(j) c * j^(-r)满足 sum(j * p(j)) 1 j np.arange(1, n1, dtypefloat) raw j ** (-r) c 1.0 / np.sum(j * raw) # 用均值约束反推归一化常数 p c * raw return p def milstein_step(S, I, R, alpha, delta, omega, p, dt, k_mean1.3418, b0.1, d0.1, sigma10.2, sigma20.1, sigma30.15): 单步 Milstein 迭代返回 (S_next, I_next, R_next) # 接触项n 层网络的加权平均感染密度 contact k_mean * np.sum(p * I) # 确定性漂移项 dS_det b - alpha * k_mean * S * contact - d * S delta * R dI_det alpha * k_mean * S * contact - (omega d) * I dR_det omega * I - (delta d) * R # 三个独立的标准正态随机数 nu1, nu2, nu3 np.random.normal(size3) sqrt_dt np.sqrt(dt) # Milstein 修正项sigma^2/2 * S * dt * (nu^2 - 1) S_next S dt * dS_det sigma1 * S * sqrt_dt * nu1 \ 0.5 * sigma1**2 * S * dt * (nu1**2 - 1) I_next I dt * dI_det sigma2 * I * sqrt_dt * nu2 \ 0.5 * sigma2**2 * I * dt * (nu2**2 - 1) R_next R dt * dR_det sigma3 * R * sqrt_dt * nu3 \ 0.5 * sigma3**2 * R * dt * (nu3**2 - 1) # 归一化投影强制 SIR1 total S_next I_next R_next S_next, I_next, R_next S_next/total, I_next/total, R_next/total # 截断到 [0,1] 避免负占比 return np.clip([S_next, I_next, R_next], 0.0, 1.0)逻辑说明compute_pj 里用均值约束 ∑j·p(j)1 反推 c这比仅做概率归一化多一个信息保证接触项的加权平均在统计意义上无偏。milstein_step 中漂移项 dS_det、dI_det、dR_det 和确定性 SIR 完全一致随机项拆成 Δt¹ᐟ² 阶和 Δt 阶两部分——前者是布朗运动的主项后者是 Milstein 独有的修正。归一化投影放在随机项之后是保证 SIR1 的关键。参数说明sigma3 资源未显式给出我取 0.15 是经验值如果模拟出的 R 序列波动比表一小适当上调dt 建议先用 0.01 试跑稳定后再尝试 0.02下面会细说。3.3 时间步长选择与噪声项注入的细节Δt 的选择直接影响稳定性。Milstein 格式的稳定性条件粗略说要求 σ²Δt 远小于 1。σ₁k0.2 时 σ²0.04Δt0.5 会产生 0.02 的二阶噪声修正看似不大但每步累积下来模拟 200 步后修正项接触总量相当于多加了 4 个单位的确定性漂移轨迹早就偏离了。我用 Δt0.01 起步每个季度对应模拟 90 步左右表一数据是按季度采样三个月约 90 天噪声项通过 np.random.normal(size3) 每步生成三个独立随机数对应三个状态变量的独立扰动。4. 参数估计从最小二乘目标函数到 Metropolis-Hastings 后验校正4.1 目标函数的设计三个 ρ² 的加权最小化给定初始值 S_{k,1}、I_{k,1}、R_{k,1}表一第一行对任意一组候选参数 (α, δ, ω)用 Milstein 格式前向推进得到每个季度的模拟值然后和表一真实值比较ρ₁² Σᵢ₌₀ⁿ (Ŝ_{k,i} - S_{k,i})²ρ₂² Σᵢ₌₀ⁿ (Î_{k,i} - I_{k,i})²ρ₃² Σᵢ₌₀ⁿ (R̂_{k,i} - R_{k,i})²总目标 J ρ₁² ρ₂² ρ₃²。这里有个值得注意的设计选择三个 ρ² 是直接相加不是加权——意味着模型假设 S、I、R 的观测误差方差相同。如果实际数据里 I 的波动明显更大通常如此可以考虑给 ρ₂² 乘权重小于 1 的系数但资源里没做属于可扩展空间。由于 Milstein 迭代里每个步都注入随机噪声同样的参数每次跑出的 J 都不同。这就是随机模型的特殊之处目标函数不是确定性的。处理办法是多次独立模拟取平均或者每个季度用同一个随机数种子做公共随机数CRN降方差。我建议先用固定种子跑主估计再换种子做敏感性分析。4.2 搜参流程与伪代码实现搜参策略是典型的非线性最小二乘从初始猜测出发观察 ρ² 变化方向向减小的方向调整 α、δ、ω。但随机噪声让 ρ² 表面粗糙不平简单梯度下降容易踩进局部极小。实际做法是多起点随机搜索加局部精修import numpy as np from scipy.optimize import minimize # 表一数据季度采样12 个时间点 data np.array([ [0.250, 0.083, 0.667], [0.250, 0.333, 0.417], [0.250, 0.250, 0.500], [0.250, 0.250, 0.500], [0.333, 0.167, 0.500], [0.333, 0.417, 0.250], [0.333, 0.250, 0.417], [0.333, 0.250, 0.417], [0.250, 0.250, 0.500], [0.333, 0.333, 0.333], [0.333, 0.417, 0.250], [0.500, 0.250, 0.250] ]) def simulate_quarterly(alpha, delta, omega, n_quarters12, steps_per_quarter90, seed42): 跑完整季度序列返回 (S_sim, I_sim, R_sim) 各 13 个点 rng np.random.default_rng(seed) p compute_pj() S, I, R data[0] # 从真实初始值出发 S_list, I_list, R_list [S], [I], [R] dt 1.0 / steps_per_quarter for _ in range(n_quarters): for _ in range(steps_per_quarter): S, I, R milstein_step(S, I, R, alpha, delta, omega, p, dt) S_list.append(S); I_list.append(I); R_list.append(R) return np.array(S_list), np.array(I_list), np.array(R_list) def objective(params, seed42): alpha, delta, omega params # 加个下限保护参数必须为正 if min(params) 0: return 1e6 S_sim, I_sim, R_sim simulate_quarterly(alpha, delta, omega, seedseed) rho1 np.sum((S_sim[1:] - data[:, 0])**2) rho2 np.sum((I_sim[1:] - data[:, 1])**2) rho3 np.sum((R_sim[1:] - data[:, 2])**2) return rho1 rho2 rho3 # 多起点搜索5 组初始猜测避免局部极小 init_guesses [(0.3, 0.2, 0.4), (0.5, 0.1, 0.3), (0.2, 0.3, 0.5), (0.4, 0.15, 0.35), (0.6, 0.05, 0.45)] best None for guess in init_guesses: res minimize(objective, guess, methodNelder-Mead, options{maxiter: 300, xatol: 1e-4, fatol: 1e-4}) if best is None or res.fun best.fun: best res print(f最优参数: alpha{best.x[0]:.3f}, delta{best.x[1]:.3f}, omega{best.x[2]:.3f}) print(f最优目标值: {best.fun:.4f})逻辑说明simulate_quarterly 从表一第一行真实值出发每季度内部跑 90 个 Milstein 步记录季度末状态objective 里把模拟值和表一后 12 行比较返回三通道残差平方和。多起点用 Nelder-Mead 是因为目标函数不光滑、无解析梯度单纯形法比梯度下降稳。每个起点独立跑 300 次迭代最后取所有起点里目标函数最小的那组。参数说明init_guesses 覆盖 α 在 0.2~0.6、δ 在 0.05~0.3、ω 在 0.3~0.5 的合理区间这个范围来自传染病 SEIR 模型的经验参数分布。如果你跑出的结果贴近边界说明数据信息不足以识别该参数需要扩大搜索范围或考虑参数固定。4.3 Metropolis-Hastings 采样给点估计补一个置信区间最小二乘只给出一组最优参数但随机模型里参数本身也是随机变量——不同参数组合可能产生几乎一样的拟合优度。这就要用摘要里提到的 Metropolis-HastingsM-H算法补一层贝叶斯推断把 ρ² 转换成一个拟似然然后从参数的后验分布采样。M-H 的口语化理解是从当前参数出发随机提出一个新参数如果新参数让目标函数更小就接受它如果更大就按概率接受——这个概率保证了采样过程最终收敛到后验分布。M-H 在传染病参数估计里的价值不是替代最小二乘而是回答α 的估计值到底有多可信。最小二乘点估计是后验的众数M-H 采样则告诉你后验分布的形状。如果后验分布很宽比如 α 的 90% 置信区间跨了 0.3说明数据对 α 的识别能力有限这时候最小二乘的最优值只是很多等价解里的一个。实现上M-H 的建议分布用对数正态扰动保证参数始终为正。每次迭代提出新参数按 Metropolis 接受率决定是否跳转。跑 20000 次采样丢弃前 5000 次作为 burn-in剩下的样本就是后验分布的近似。def mh_sampler(n_iter20000, burn_in5000, tune0.2, seed7): M-H 采样返回后验样本数组 (n_iter-burn_in, 3) rng np.random.default_rng(seed) # 以最小二乘最优解作为起点加快收敛 current np.array([0.5, 0.15, 0.4]) current_obj objective(current) samples [] for i in range(n_iter): proposal current * np.exp(tune * rng.normal(size3)) # 乘性扰动保证正数 proposal_obj objective(proposal) # 拟似然比随机模型下目标函数有噪声这里做 3 次平均平滑 if proposal_obj current_obj or rng.random() np.exp(current_obj - proposal_obj): current, current_obj proposal, proposal_obj if i burn_in: samples.append(current.copy()) return np.array(samples)逻辑说明乘性对数正态扰动天然保证提议参数非负符合 α、δ、ω 的物理约束。接受概率用 exp(ΔJ) 形式和标准 M-H 用似然比的写法等价——因为这里最小二乘目标 J 与负对数似然成正比。随机模型目标函数本身有波动我建议每次评估目标函数时固定随机种子否则接受率会被噪声主导。5. 避坑与常见问题复现这套模型时最容易翻车的 5 个点5.1 状态变量出现负值或超过 1现象模拟中途 S 或 I 变成负数或者 R 超过 1后续迭代直接发散。原因是 Milstein 修正项 (σ²/2)SΔt(ν²-1) 里的 ν² 可以很大当 Δt 不够小时修正项幅度超过漂移项把状态变量推出物理边界。解决每步做完随机修正后用 np.clip(…, 0, 1) 截断再归一化。这不算漂亮但很实用截断引入的偏差在 Δt→0 时消失。5.2 目标函数波动导致搜参不收敛现象Nelder-Mead 迭代到后期目标函数在两个相近的参数点之间来回跳不下降。原因是每步模拟都重新生成随机数目标函数本身是随机变量我见过的最惨案例是同一组参数两次运行目标值差 5 倍。解决在 objective 函数里固定 seed或者每个参数做 3~5 次独立模拟取平均目标值。固定 seed 会损失部分随机性但搜参阶段要的是目标函数可比较。5.3 参数不可识别多组参数值拟合效果几乎一样现象两个初始猜测收敛到完全不同的参数但目标函数值接近。原因SIR 模型里 α、δ、ω 存在补偿效应——α 增大、ω 减小可能产生相似的感染峰值。我用 M-H 采样看过后验分布α 和 ω 在二维平面上的等高线是拉长的椭圆说明数据信息不足以单独识别这两个参数。解决固定其中一个比如 δ0.1专门扫另外两个或者在报告结果时给出后验分布而不仅是点估计。5.4 p(j) 归一化条件用错现象接触项 Σjp(j)I_j 算出的均值明显偏大模拟的感染人数系统性偏高。原因只做了概率归一化 Σp(j)1没有做均值约束 Σjp(j)1。用资源里的 r3、n5两种归一化给出的 c 值差约 30%足以改变模型动态。解决compute_pj 函数里用均值约束求 c每次改 n 或 r 后先打印验证 Σjp(j)≈1。5.5 季度内步数太少导致离散化偏差混入参数估计现象用 Δt1一个季度一步跑出来的最优参数明显偏离合理范围。原因Milstein 格式在 Δt 过大时失去收敛性离散化偏差直接进了参数反演。解决先固定 Δt0.01 跑通流程再试 0.005 和 0.02看参数估计结果是否稳定——如果 α 在两个 Δt 下差超过 20%说明离散化还没收敛要继续缩小步长。6. 收尾技巧模型自洽性检验与后验诊断一套参数估计做出来最怕的不是不收敛而是收敛到一个自洽性欠佳的模型——模拟分布和真实数据差异大但目标函数碰巧很小。我习惯用三个验证步骤给结果兜底。第一步是验证总量守恒。跑完整模拟后输出 SIR 的残差序列理论应为零。如果残差在某个区间集中偏离比如总是 0.02说明归一化投影和截断之间的顺序有问题需要回看代码里投影是否放在截断之后。第二步是比较模拟分布与真实数据的统计量。对最优参数跑 100 次独立模拟得到每个季度 S、I、R 的均值 ± 标准差带检查真实数据是否落在这个带内。占比类传染病数据80% 以上的点落在一个标准差带内是比较合理的水平。如果真实值频繁超出两个标准差带说明 σ 参数设置偏小模型不确定性不足以覆盖观测波动要上调 σ。第三步是参数后验诊断。用 M-H 采样结果计算每个参数的 Gelman-Rubin 统计量多链收敛指标或者直接看后验直方图的形状系数。如果后验分布接近高斯说明数据信息充分如果后验分布出现双峰说明存在两组参数同样解释数据需要回到实际业务场景里判断哪一组更合理。我在复现这套资源时印象最深的一个教训是最小二乘给出的最优参数看着挺好但 M-H 后验分布却宽得吓人——这说明数据对参数的约束能力远低于直觉判断。从那以后每次做随机传染病模型的参数估计我都强制走一遍固定种子搜参 → 换种子验证稳定性 → M-H 采样看后验宽度的流程先做完这套再谈结论。希望帮到你。本文还有配套的精品资源点击获取