简介这份资源围绕传染病动力学建模中的随机微分方程展开面向流行病建模、公共卫生干预评估方向的研究人员及高校研究生帮助读者掌握环境噪声下疫情传播的随机建模思路与数据驱动的参数估计方法。压缩包共1个文件为docx文档整体约8.09MB内容涵盖模型推导、离散化数值格式与参数估计流程。文档以标准布朗运动刻画不可预见变化构建可描述真实疫情动态的随机传染模型并基于2021年初至2022年底共21个月的中国官方统计数据涉及确诊病例、治愈人数与疫苗接种人次等维度。参数估计部分采用Metropolis-Hastings算法完成贝叶斯推断下的最大似然估计验证大样本抽样解的稳定性并给出Milstein格式等数值仿真关键细节附有离散化方程、参数取值与季度数据表。已有151人学习适合作为随机微积分与传染病建模交叉研究的实操参考。1. 传染病动力学建模中的随机微分方程从确定性SIR到带噪声的参数估计做过传染病建模的人大多有过这种体验用确定性 SIR 模型拟合早期数据曲线看着挺漂亮可一旦把预测窗口拉长置信区间窄得离谱现实数据却总在区间外游走。问题不在参数拟合得不好而在于确定性模型把传播过程当成了钟表——给定初值和参数未来完全确定。真实传播里接触率随人群行为波动、超级传播事件随机出现、检测和报告存在延迟噪声这些随机性不是误差而是系统本身的一部分。随机微分方程SDE建模就是把这份随机性写进方程让参数估计从「点估计」变成「分布估计」。这篇笔记面向已经会跑确定性模型、想进一步做随机版本参数估计的从业者从 SDE 的构造、数值求解、到基于实际疫情数据的参数反演把能复现的路径和踩过的坑一次讲清。2. 把随机性写进传播方程从 SIR 到 Itô 形式的三种构造方式2.1 为什么确定性 SIR 在参数估计上会「过度自信」确定性 SIR 的方程是 dS/dt -βSI/NdI/dt βSI/N - γI。它的解是一条光滑曲线拟合时最小化残差平方和得到的 β 和 γ 是使曲线最贴近观测点的值。但观测数据本身是随机过程的一次实现不同随机种子下同一组参数会产生不同轨迹。确定性模型没有机制表达这种「同一参数、不同轨迹」的离散程度于是拟合出的参数置信区间只反映观测误差不反映过程噪声。结果就是预测区间偏窄参数估计的方差被系统性低估。常见做法是在确定性方程上直接加高斯噪声写成 dX f(X)dt σ dW。这看起来简单但有个隐蔽问题加性噪声不随状态变化当感染人数接近零时噪声可能把 I 推成负数。更合理的是乘性噪声让波动幅度与当前状态成比例。2.2 三种构造加性噪声、乘性噪声与跳扩散实际建模中我一般按数据特征选构造方式。早期数据稀疏、报告延迟明显时加性噪声够用感染规模较大、行为波动显著时乘性噪声更贴合出现超级传播事件或聚集性爆发时跳扩散过程能捕捉那些「突然跳一下」的轨迹。乘性噪声的 Itô 形式写成import numpy as np def sir_sde_multiplicative(state, t, beta, gamma, sigma): SIR 模型的乘性噪声 SDE state: [S, I, R] sigma: 噪声强度控制波动幅度 返回漂移项和扩散项 S, I, R state N S I R # 漂移项确定性部分 drift_S -beta * S * I / N drift_I beta * S * I / N - gamma * I drift_R gamma * I # 扩散项乘性噪声波动与 sqrt(I) 成比例 diff_S -sigma * S * I / N diff_I sigma * S * I / N diff_R sigma * I return np.array([drift_S, drift_I, drift_R]), np.array([diff_S, diff_I, diff_R])这段代码里漂移项就是确定性 SIR 的右端扩散项则把噪声强度 σ 乘上状态相关的量。σ 不是随便设的它决定了轨迹的离散程度。σ 太大轨迹会频繁穿越零附近数值求解不稳定σ 太小随机性体现不出来和确定性模型没区别。我一般先用 σ 在 0.01 到 0.1 之间试看模拟轨迹的包络是否覆盖观测数据的波动范围。跳扩散的构造更复杂需要额外定义跳跃强度和跳跃幅度分布。常见做法是泊松跳跃强度 λ 控制单位时间跳跃次数幅度用对数正态或指数分布。跳扩散适合有明确聚集性事件的数据但参数多、估计难新手建议先从乘性噪声入手。2.3 数值求解Euler-Maruyama 与 Milstein 的取舍SDE 没有解析解必须数值求解。最常用的是 Euler-Maruyama 方法它对应常微分方程里的 Euler 法def euler_maruyama(state0, t_span, drift_fn, diff_fn, params, dt0.01): Euler-Maruyama 求解 SDE state0: 初值 t_span: (t_start, t_end) drift_fn, diff_fn: 漂移和扩散函数 params: 模型参数 dt: 时间步长 t_start, t_end t_span n_steps int((t_end - t_start) / dt) states np.zeros((n_steps 1, len(state0))) states[0] state0 for i in range(n_steps): t t_start i * dt drift, diff drift_fn(states[i], t, *params) dW np.random.normal(0, np.sqrt(dt), sizelen(state0)) states[i1] states[i] drift * dt diff * dW return statesdt 的选择很关键。dt 太大离散化误差会累积轨迹偏离真实 SDE 的解dt 太小计算量爆炸。我一般先用 dt0.01 跑一遍再把 dt 减半看结果是否稳定。如果两次结果差异明显说明 dt 还不够小。Milstein 方法在扩散项依赖状态时精度更高但需要计算扩散项的导数实现复杂。对于 SIR 这类扩散项相对简单的模型Euler-Maruyama 在 dt0.001 到 0.01 之间通常够用。提示Euler-Maruyama 的收敛阶是 0.5强收敛比确定性 Euler 法的 1 阶低。这意味着要达到同样精度SDE 需要更小的步长或更多的样本路径。3. 基于实际疫情数据的参数估计从最小二乘到贝叶斯反演3.1 数据预处理报告延迟与检测窗口的修正实际疫情数据不是干净的 I(t)而是每日新增报告病例。报告延迟意味着今天确诊的人可能几天前就感染了检测窗口则决定了有多少感染者被漏掉。直接拿报告数当 I(t) 拟合参数估计会有系统性偏差。我一般做两步修正。第一步用滑动平均或样条插值把离散的每日报告转成连续曲线减少报告日的周期性波动。第二步假设报告延迟服从某个分布常用伽马分布或对数正态用反卷积把报告数还原成感染发生数。这一步不需要很精确但能显著改善 β 的估计。import pandas as pd from scipy.ndimage import gaussian_filter1d def preprocess_cases(case_series, sigma2.0): 对每日报告病例做平滑处理 case_series: 每日新增报告数 sigma: 高斯平滑核宽度 smoothed gaussian_filter1d(case_series.astype(float), sigmasigma) # 保证非负 smoothed np.maximum(smoothed, 0) return smoothedsigma 控制平滑程度。sigma 太小噪声没压住sigma 太大早期增长信号被抹平。我一般从 sigma1 开始逐步加到 3看拟合残差是否还有明显结构。如果残差图呈现周期性波动说明平滑不够如果峰值被削平说明平滑过度。3.2 最小二乘估计目标函数与优化器选择有了平滑后的感染曲线下一步是找一组参数 (β, γ, σ) 使模型模拟的轨迹最贴近数据。最小二乘的目标函数是模拟轨迹与观测轨迹的加权平方误差。由于 SDE 每次模拟结果不同目标函数本身是随机的不能直接用梯度下降。常见做法是固定随机种子对同一组参数多次模拟取平均再用优化器搜索。优化器我一般用 Nelder-Mead 或差分进化它们不需要梯度对随机目标函数更鲁棒。from scipy.optimize import differential_evolution def loss_function(params, observed_I, state0, t_span, n_sim10): 计算参数对应的损失 params: [beta, gamma, sigma] observed_I: 观测感染曲线 n_sim: 每个参数组合模拟次数 beta, gamma, sigma params losses [] for seed in range(n_sim): np.random.seed(seed) states euler_maruyama(state0, t_span, sir_sde_multiplicative, None, (beta, gamma, sigma), dt0.01) simulated_I states[:, 1] # 只比较观测窗口内的部分 n_obs len(observed_I) loss np.mean((simulated_I[:n_obs] - observed_I) ** 2) losses.append(loss) return np.mean(losses) # 参数搜索范围 bounds [(0.1, 1.0), (0.05, 0.5), (0.001, 0.2)] result differential_evolution(loss_function, bounds, args(observed_I, state0, (0, len(observed_I))), maxiter100, tol1e-6)n_sim 的选择是个权衡。n_sim 太小损失函数噪声大优化器容易卡在局部n_sim 太大计算时间成倍增加。我一般先用 n_sim5 粗搜找到大致范围后再用 n_sim20 精搜。differential_evolution 的 maxiter 和 tol 控制收敛maxiter 太小结果不稳定太大浪费时间。tol 设 1e-6 通常够用。3.3 贝叶斯参数估计先验选择与后验采样最小二乘给的是点估计贝叶斯方法给的是后验分布。对于 SDE 模型似然函数没有解析形式必须用近似贝叶斯计算ABC或粒子滤波。ABC 的思路是从先验采样参数模拟轨迹如果模拟轨迹与观测的距离小于阈值 ε就接受这个参数样本。def abc_rejection(observed_I, state0, t_span, prior_ranges, n_samples10000, epsilon0.5): ABC 拒绝采样 prior_ranges: 每个参数的先验范围 epsilon: 接受阈值 accepted [] for _ in range(n_samples): beta np.random.uniform(*prior_ranges[beta]) gamma np.random.uniform(*prior_ranges[gamma]) sigma np.random.uniform(*prior_ranges[sigma]) np.random.seed(np.random.randint(0, 10000)) states euler_maruyama(state0, t_span, sir_sde_multiplicative, None, (beta, gamma, sigma), dt0.01) simulated_I states[:len(observed_I), 1] distance np.sqrt(np.mean((simulated_I - observed_I) ** 2)) if distance epsilon: accepted.append((beta, gamma, sigma)) return np.array(accepted)epsilon 是 ABC 的核心参数。epsilon 太大接受率高但后验近似粗糙epsilon 太小接受率极低计算量爆炸。我一般先跑一小批样本看距离的分布取距离的 5% 分位数作为 epsilon。先验范围的选择也很关键β 的先验通常取 0.1 到 1.0γ 取 0.05 到 0.5σ 取 0.001 到 0.2。先验太宽接受率低先验太窄可能把真值排除在外。注意ABC 拒绝采样在高维参数空间效率很低。如果参数超过 3 个建议改用 ABC-SMC 或粒子滤波否则计算时间不可接受。4. 避坑与排查SDE 参数估计里最容易翻车的五个地方4.1 轨迹穿越零导致数值爆炸现象模拟过程中 I 变成负数后续轨迹完全失控损失函数返回 NaN。原因乘性噪声在 I 接近零时扩散项也接近零但 Euler-Maruyama 的离散化误差可能让 I 一步跨到负值。一旦 I 为负漂移项和扩散项都失去物理意义。解决在每步更新后加一个截断把 I 限制在非负范围。更优雅的做法是用对数变换令 x log(I)在 log 空间做 SDE 模拟保证 I 始终为正。但对数变换会改变噪声结构需要重新推导扩散项。4.2 噪声强度 σ 与观测误差混淆现象拟合出的 σ 很大轨迹波动剧烈但预测区间仍然不覆盖未来数据。原因σ 同时吸收了过程噪声和观测误差。如果观测数据本身噪声大σ 会被高估导致模型把观测误差当成过程噪声预测时过度离散。解决如果有多条独立数据流如不同地区的报告可以联合拟合把共享的 σ 和地区特定的观测误差分开。单条数据流时只能通过先验约束 σ 的范围避免它无限增大。4.3 优化器陷入局部最优现象不同随机种子下 differential_evolution 给出的参数差异很大损失函数值也差一个量级。原因SDE 的损失函数是随机的且可能存在多个局部极小。参数空间里 β 和 γ 存在相关性沿某些方向损失变化平缓优化器容易停在平坦区域。解决增加 n_sim降低损失函数噪声扩大参数搜索范围避免边界卡住用多组初值分别优化取损失最小的那组。我一般会跑 3 到 5 次独立优化看结果是否收敛到同一区域。4.4 时间步长 dt 与观测间隔不匹配现象模拟轨迹在观测点之间的行为与数据不符拟合残差呈现系统性偏差。原因dt 太小模拟轨迹的分辨率远高于观测数据损失函数只比较观测点忽略了观测间隔内的动态。dt 太大离散化误差累积轨迹偏离真实 SDE 解。解决让 dt 至少比观测间隔小一个量级。如果数据是每日报告dt 取 0.01 天约 15 分钟通常够用。同时检查模拟轨迹在观测点之间的插值是否合理必要时对模拟轨迹做与观测相同的平滑处理。4.5 先验信息与数据冲突现象贝叶斯后验分布集中在先验边界或者后验与先验几乎重合数据没有提供有效信息。原因先验范围设置不当把真值排除在外或者数据信息量不足似然函数太平坦后验被先验主导。解决先用最小二乘或 ABC 粗搜确定参数的大致范围再据此设置先验。如果后验集中在边界说明先验范围需要放宽。如果后验与先验重合说明数据对该参数不敏感需要考虑增加数据量或简化模型。5. 进阶技巧用粒子滤波做在线参数估计与模型选择5.1 粒子滤波为什么适合 SDE 参数估计SDE 的似然函数没有解析形式但粒子滤波可以通过序贯蒙特卡洛近似似然。它的核心思想是用一组带权重的粒子表示状态分布每步根据观测更新权重重采样避免退化。对于参数估计可以把参数也当成状态的一部分做联合估计。粒子滤波的优势在于在线更新新数据到来时不需要重新拟合全部历史天然处理非线性非高斯可以同时估计状态和参数。缺点是计算量大粒子数不够时权重退化严重。def particle_filter(observations, n_particles1000, param_boundsNone): 粒子滤波做 SDE 参数在线估计 observations: 观测序列 n_particles: 粒子数 param_bounds: 参数范围 # 初始化粒子状态和参数联合采样 particles np.zeros((n_particles, 5)) # S, I, R, beta, gamma weights np.ones(n_particles) / n_particles for i in range(n_particles): particles[i, 0] 990 # S0 particles[i, 1] 10 # I0 particles[i, 2] 0 # R0 particles[i, 3] np.random.uniform(*param_bounds[beta]) particles[i, 4] np.random.uniform(*param_bounds[gamma]) estimates [] for t, obs in enumerate(observations): # 预测步用 SDE 推进每个粒子 for i in range(n_particles): S, I, R, beta, gamma particles[i] N S I R dW np.random.normal(0, np.sqrt(0.01)) sigma 0.05 S_new S - beta * S * I / N * 0.01 - sigma * S * I / N * dW I_new I (beta * S * I / N - gamma * I) * 0.01 sigma * S * I / N * dW R_new R gamma * I * 0.01 sigma * I * dW particles[i] [max(S_new, 0), max(I_new, 0), max(R_new, 0), beta, gamma] # 更新步根据观测计算权重 for i in range(n_particles): pred_I particles[i, 1] weights[i] * np.exp(-0.5 * ((obs - pred_I) / 10) ** 2) weights / np.sum(weights) # 重采样 if 1.0 / np.sum(weights ** 2) n_particles / 2: indices np.random.choice(n_particles, n_particles, pweights) particles particles[indices] weights np.ones(n_particles) / n_particles # 记录当前估计 beta_est np.sum(weights * particles[:, 3]) gamma_est np.sum(weights * particles[:, 4]) estimates.append((beta_est, gamma_est)) return estimates这段代码里粒子数 n_particles 决定估计精度。粒子太少权重退化快估计不稳定粒子太多计算慢。我一般从 1000 开始看有效粒子数1/sum(weights^2)是否经常低于 n_particles/2如果是就加倍。观测噪声的标准差这里设 10需要根据数据调整设得太小会导致权重极端集中设得太大则观测信息利用不足。5.2 模型选择加性噪声还是乘性噪声有了参数估计方法下一个问题是选哪个模型。加性噪声和乘性噪声给出不同的后验分布哪个更可信我一般用信息准则或贝叶斯因子比较。信息准则方面WAIC广泛可用信息准则适合贝叶斯模型它同时考虑拟合优度和有效参数数。计算 WAIC 需要后验样本的对数似然粒子滤波可以输出每步的似然估计。贝叶斯因子则需要计算边际似然粒子滤波的似然估计可以累乘得到。实际操作中我一般先跑两种噪声结构的参数估计比较后验预测区间对观测数据的覆盖率和区间宽度。乘性噪声如果显著改善覆盖率且区间没有过度膨胀就选乘性如果两者差不多选加性因为参数少、计算快。5.3 一个具体技巧用轮廓似然确定参数可辨识性SDE 参数估计里β 和 γ 往往存在相关性单独看后验边缘分布可能都正常但联合分布呈香蕉形。轮廓似然能揭示这种相关性固定 β 在一系列值上对每个 β 优化 γ 和 σ画出损失函数随 β 的变化曲线。如果曲线平坦说明 β 不可辨识如果曲线有明确极小说明 β 可辨识。def profile_likelihood(beta_values, observed_I, state0, t_span): 轮廓似然固定 beta优化 gamma 和 sigma profile [] for beta in beta_values: def loss(params): gamma, sigma params return loss_function([beta, gamma, sigma], observed_I, state0, t_span, n_sim5) res differential_evolution(loss, [(0.05, 0.5), (0.001, 0.2)], maxiter50) profile.append(res.fun) return np.array(profile)beta_values 的范围根据先验或粗搜结果确定一般取 20 到 50 个点。每个点上的优化用 differential_evolution 可能太慢可以改用 Nelder-Mead因为只有两个参数。轮廓似然曲线如果出现多个极小说明模型可能存在多解需要检查数据是否覆盖了足够多的动态阶段。我自己的习惯是每次做 SDE 参数估计先跑轮廓似然看参数可辨识性再决定是否值得做完整的贝叶斯推断。如果轮廓似然平坦贝叶斯后验也会很宽这时候与其花时间采样不如回头检查数据质量或简化模型。这个习惯帮我省下了不少算力也避免了很多「后验不收敛」的玄学问题。希望帮到你。本文还有配套的精品资源点击获取