SDE传染病建模实战：新冠数据参数估计与避坑指南

简介传染病动力学建模中随机微分方程能把环境扰动引入传播过程比确定性模型更贴近真实疫情波动。资源围绕随机SIR型传染病模型展开面向流行病建模研究人员、公共卫生政策制定者及高校研究生旨在解决随机建模、数据驱动参数估计与数值仿真实现等问题。资源包仅含1个DOCX文档约8.09MB内容涵盖随机微分方程建模、布朗运动驱动项设定、疫情数据参数估计流程以及M-H采样与Milstein格式离散化的关键推导并给出参数取值、递推公式和S/I/R状态比例表。文档进一步展示了离散化过程、噪声强度设置及参数自适应调整细节便于复现M-H抽样下的贝叶斯参数推断也可作为教学辅助材料。目前已有151人浏览学习适合需要将随机微积分工具应用于传染病数据分析的研究者作为进阶参考。1. 随机微分方程进入传染病建模新冠疫情数据逼出来的三个选择传染病动力学建模如果只靠确定性 SEIR 拟合新冠疫情每日新增数据第一个星期就会撞墙模型输出是光滑曲线真实数据是剧烈抖动的锯齿R0 今天算出来 3.2下周同一批数据重算就掉到 2.4。墙后面是随机微分方程SDE——在 ODE 右侧加一个扩散项让趋势和涨落同时被建模。这篇笔记只讲两件事SDE 模型怎么构造以及用新冠疫情数据做参数估计哪条路线最稳适合已跑通 SEIR、正被预测区间过窄和参数漂移折磨的建模工程师。我把三个最关键的决策点放在前面噪声加在哪个状态变量上、用 Itô 还是 Stratonovich 积分、观测模型用高斯还是负二项。这三件事定下来后面的参数估计方法才能稳定输出可用结果。代码和参数表部分可以直接照着改坑的部分是我真实踩过的不是教科书推演。2. 传染病 SDE 模型的方程构造噪声项加在哪、用哪种积分2.1 从 SEIR 到 SDE先分清“人口学噪声”和“传播噪声”标准 SEIR 的 ODE 形式是dS/dt −βSI/NdE/dt βSI/N − σEdI/dt σE − γIdR/dt γI这里 S、E、I、R 是易感、潜伏、感染、恢复人数N 是总人口。确定性模型把所有个体当成连续的“平均人”每个时刻的转移速率是确定的。但真实传播不是这样一个人一天接触多少人、接触后是否真的被传染都是随机事件。当感染人数只有几百人时这些随机性的相对影响非常大甚至可能让疫情早期自然灭绝。常见的 SDE 化方式有两类。第一类是人口学噪声化学 Langevin 近似。把 S→E、E→I、I→R 三个转移事件看成速率分别为 λ₁βSI/N、λ₂σE、λ₃γI 的马尔可夫跳过程当群体规模足够大时扩散近似给出dS −λ₁dt √λ₁ dW₁dE (λ₁ − λ₂)dt √λ₁ dW₁ − √λ₂ dW₂dI (λ₂ − λ₃)dt √λ₂ dW₂ − √λ₃ dW₃dR λ₃dt √λ₃ dW₃注意 S 和 E 的噪声项共享同一个 dW₁因为一个人离开 S 必然进入 E相关性是守恒出来的。第二类是传播噪声。把有效接触率 β 看作围绕均值波动的随机过程扩散项直接乘在感染压力上dS −βSI/N dt − σ_β(SI/N)dW传播噪声描述的是环境因素温度、湿度、人群聚集政策的松紧带来的整体扰动所有个体同时受影响。我个人的工程习惯是用传播噪声模型做参数拟合因为它的参数少、数值稳定拟合完再跑一个人口学噪声模型做灭绝概率和早期风险的敏感性检查两个模型的噪声项含义不同不能混着用。2.2 Itô 还是 Stratonovich传染病模型里为什么默认选 Itô同一个 SDE用 Itô 积分和 Stratonovich 积分会得到不同的样本路径和不同的参数含义。区别在于积分时取区间左端点还是中点。对传染病模型这个选择不是纯数学偏好它影响“β 到底代表什么”。如果噪声来自过程本身人口学噪声那么连续时间极限下自然得到 Itô 积分没有选择余地。如果噪声来自外界环境比如温度波动影响传播率Stratonovich 在物理上更自然因为它对应的白噪声是随机环境的连续近似。但从离散数据做参数估计时我们用 Euler-Maruyama 离散化它天然是 Itô 积分所以实践中几乎全部论文和代码都默认 Itô。有一条需要注意Stratonovich 解释下的漂移项等于 Itô 漂移加上一个修正项 (1/2)g(x)g′(x)。对于传播噪声模型 g(S)σ_β·βSI/N这个修正项的量级是 (1/2)(σ_β·βI/N)²在疫情早期 I/N 很小的时候可以忽略但在感染比例高、噪声强度大的时候会把 R0 的估计往高推几个百分点。我的做法是固定用 Itô但在报告里写明如果审稿人或同事问起来把修正项算一遍确认影响小于预期误差就再解释。2.3 三类常见 SDE 模型的适用场景与参数表在真实项目里我见过并且用过三类噪声结构各有各的适用边界。噪声模型扩散项形式适用场景主要坑化学 Langevin√λ₁ √λ₂ √λ₃小群体、早期爆发、灭绝概率大群体下步长必须很小计算慢传播率乘性噪声σ_β·βSI/N拟合现实疫情曲线的波动σ_β 和 β 强相关可辨识性差参数随机游走dβκ(β̄−β)dtσ_βdW政策干预导致 β 时变参数多短观测窗口易过拟合模型参数表是每次估计前要确定的我不建议让优化器一口气全估计出来。参数含义典型范围来源β有效接触率/天0.1—1.5估计σ潜伏期倒数 E→I/天1/3—1/7文献固定γ恢复率 I→R/天1/5—1/10文献固定σ_β传播噪声强度0.01—0.3估计ρ病例报告率0.1—1.0独立血清学调查φ观测离散参数负二项3—20估计或固定下面是一个能直接跑的 SEIR-SDE 模拟器Euler-Maruyama 积分支持两种噪声模型。代码里我把状态截断在 0避免噪声把人数推成负数。import numpy as np def seir_sde_simulate(beta, sigma_e, gamma, S0, E0, I0, R0, noise_typetransmission, sigma_noise0.1, t_end120, dt0.05, seed0): Euler-Maruyama 积分 SEIR 随机微分方程。 参数说明 ---- beta : 有效接触率/天 sigma_e : E→I 速率/天 gamma : I→R 速率/天 noise_type : transmission 为传播噪声g sigma_noise * beta*S*I/N langevin 为化学 Langevin 噪声g sqrt(rate) sigma_noise : 传播噪声强度仅 transmission 模式使用 rng np.random.default_rng(seed) n_steps int(t_end / dt) t np.linspace(0, t_end, n_steps 1) S np.empty(n_steps 1) E np.empty(n_steps 1) I np.empty(n_steps 1) R np.empty(n_steps 1) S[0], E[0], I[0], R[0] S0, E0, I0, R0 N S0 E0 I0 R0 for k in range(n_steps): s, e, i S[k], E[k], I[k] lam1 beta * s * i / N if N 0 else 0.0 lam2 sigma_e * e lam3 gamma * i if noise_type langevin: dW1 rng.normal(0, np.sqrt(dt)) dW2 rng.normal(0, np.sqrt(dt)) dW3 rng.normal(0, np.sqrt(dt)) ds -lam1 * dt np.sqrt(max(lam1, 0.0)) * dW1 de (lam1 - lam2) * dt np.sqrt(max(lam1, 0.0)) * dW1 \ - np.sqrt(max(lam2, 0.0)) * dW2 di (lam2 - lam3) * dt np.sqrt(max(lam2, 0.0)) * dW2 \ - np.sqrt(max(lam3, 0.0)) * dW3 dr lam3 * dt np.sqrt(max(lam3, 0.0)) * dW3 else: dW rng.normal(0, np.sqrt(dt)) g sigma_noise * lam1 ds -lam1 * dt - g * dW de (lam1 - lam2) * dt g * dW di (lam2 - lam3) * dt dr lam3 * dt S[k 1] max(S[k] ds, 0.0) E[k 1] max(E[k] de, 0.0) I[k 1] max(I[k] di, 0.0) R[k 1] max(R[k] dr, 0.0) return t, S, E, I, R逻辑说明每一小步都把四个状态按转移速率推进传播噪声模式下 dW 只加在 S 和 E 上因为 I 和 R 的转移速率本身没有环境噪声化学 Langevin 模式则严格保持 S→E→I→R 的守恒关系噪声用 √rate 缩放保证大速率事件有大的随机涨落。步长 dt 设为 0.05 天对应每天 20 个中间子步这是传播噪声模型下精度和速度的折中。参数说明sigma_noise 的典型量级是 0.05—0.2再大容易让 S 触底seed 参数必须在敏感性分析时固定否则无法区分随机波动和参数变化。如果想换 I 与 R 的噪声结构把 di、dr 两行改成与传播噪声共享 dW 即可但那样会让峰值预测飘得很难看。3. 用新冠疫情数据做参数估计伪极大似然与贝叶斯 MCMC3.1 离散观测下的似然函数用 Euler-Maruyama 搭桥SDE 的转移密度几乎没有解析形式这是参数估计最大的障碍。好消息是我们用 Euler-Maruyama 做仿真它也给出了转移密度的近似——给定 X_t下一时刻 X_{tΔt} 近似服从高斯分布均值是 X_t f(X_t)Δt方差是 g(X_t)²Δt多维情形是 GGᵀΔt。把观测时间点之间的多步转移乘起来就得到一个伪似然函数。注意“伪”字只有当 Δt 趋近于 0 时这个高斯逼近才严格成立。对每日疫情数据来说观测间隔是 1 天直接用 1 天做单步高斯逼近误差太大。所以我在构造似然时总是把 1 天拆成 20 个 0.05 天的内部子步只在观测时刻与数据比较。这样相当于把数值离散误差往下压了一到两个量级。还有一个工程问题每日新增病例是流 S→E 的积分不是任何单独一个状态变量。因此输出端要取 S 的逐日变化量 −(S_{t1}−S_t)而不是直接拿 I 和观测值比。报告延迟用参数 delay_days 对齐第 t 天报告的病例对应第 t−delay_days 天的感染事件。3.2 伪极大似然估计代码与参数设置下面这段代码把 SEIR-SDE 模拟器封装成每日新增感染序列再用负二项观测模型构造伪对数似然交给 scipy 的 Nelder-Mead 做最大化。负二项能处理计数数据的过离散比高斯假设稳健得多。from scipy.optimize import minimize from scipy.stats import nbinom import numpy as np def seir_sde_daily_incidence(beta, gamma, y_obs, S0, I0, population, sigma_e1/3, delay_days6, dt0.05, sigma_noise0.1, seed0): 模拟 SEIR-SDE返回对齐报告延迟后的每日新增感染。 t_end len(y_obs) delay_days 2 t, S, E, I, R seir_sde_simulate( beta, sigma_e, gamma, S0, I0, 0.0, 0.0, noise_typetransmission, sigma_noisesigma_noise, t_endt_end, dtdt, seedseed) # 按天抽样计算每日 S 的减少量作为新感染 S_daily S[::int(1 / dt)] daily_inf np.maximum(S_daily[:-1] - S_daily[1:], 0) aligned daily_inf[delay_days:delay_days len(y_obs)] return aligned def neg_loglik_beta_gamma(params, y_obs, S0, I0, population, rho0.3, phi5.0, delay_days6): beta, gamma params if beta 0 or gamma 0: return 1e10 y_sim seir_sde_daily_incidence(beta, gamma, y_obs, S0, I0, population, delay_daysdelay_days) mu np.maximum(rho * y_sim, 1e-6) # 负二项均值 mu离散参数 phi方差 mu mu^2/phi p phi / (phi mu) return -np.sum(nbinom.logpmf(np.round(y_obs).astype(int), phi, p)) # 假设 y_obs 是清洗后的每日新增病例序列 y_obs np.array([5, 8, 12, 15, 23, 28, 36, 51, 62, 80, ...]) # 实际数据 population 8_000_000 S0 population - 1 I0 10.0 res minimize( lambda params: neg_loglik_beta_gamma(params, y_obs, S0, I0, population), x0[0.8, 0.14], methodNelder-Mead, options{maxiter: 2000, xatol: 1e-4, fatol: 1e-3}) print(beta , res.x[0], gamma , res.x[1])逻辑说明minimize 接受的是一维参数向量 [beta, gamma]负对数似然越小越好。每天内部跑 20 个 Euler-Maruyama 子步因此 optimize 每轮迭代都要完整模拟一遍疫情过程参数接近真实值时拟合优度迅速上升形成平滑的优化面。Nelder-Mead 不需要梯度对 SDE 模拟这种带随机性的目标函数很稳换 L-BFGS-B 的话随机性会让梯度爆炸。参数说明rho 是报告率这个例子固定成 0.3。phi 是负二项离群参数固定 5 表示观测方差大约是均值的 1.2 倍。如果数据噪声特别大把 phi 设小一点但注意 phi 与 sigma_noise 有部分重叠两个都放开容易导致不可辨识。I0 我建议先取早期平均日增病例除以 γρ算出来再代入不要直接从 1 开始扫。3.3 贝叶斯 MCMC先验、Metropolis 采样与收敛诊断伪极大似然给的是点估计但疫情决策需要区间。把伪似然当成观测模型配上参数先验就能用 MCMC 采出后验分布。下面给出一个手写随机游走 Metropolis 的完整例子用来估计 beta 和 gamma。这段代码刻意不引第三方贝叶斯库因为 SDE 的伪似然在 PyMC 里写自定义函数时绕不开直接手写反而能看清每一步在做什么。def log_posterior(params, y_obs, S0, I0, population, prior_beta(0.5, 0.5), prior_gamma(0.2, 0.2)): beta, gamma params if not (0.0 beta 3.0 and 0.0 gamma 1.0): return -np.inf nll neg_loglik_beta_gamma(params, y_obs, S0, I0, population) # 对数正态先验 lp (-0.5 * ((np.log(beta) - np.log(prior_beta[0])) / prior_beta[1])**2 - 0.5 * ((np.log(gamma) - np.log(prior_gamma[0])) / prior_gamma[1])**2) return -nll lp def metropolis(log_posterior, init, n_iter20000, step0.03, seed1): rng np.random.default_rng(seed) cur np.array(init, dtypefloat) cur_lp log_posterior(cur) samples np.zeros((n_iter, 2)) accept 0 for i in range(n_iter): prop cur rng.normal(0, step, size2) prop_lp log_posterior(prop) if np.log(rng.uniform()) prop_lp - cur_lp: cur, cur_lp prop, prop_lp accept 1 samples[i] cur return samples, accept / n_iter samples, acc_rate metropolis( lambda p: log_posterior(p, y_obs, S0, I0, population), init[0.8, 0.14], n_iter20000, step0.03) print(acceptance rate:, acc_rate) print(posterior mean beta:, samples[:, 0].mean()) print(posterior mean gamma:, samples[:, 1].mean())逻辑说明MH 每一步从当前参数附近走一步接受概率由后验比决定步长 step 控制探索半径。20%—40% 的接受率是随机游走 MH 的经验区间太高说明步子太小太低说明步子太大。参数说明先验我把 beta 的中心放在 0.5、gamma 放在 0.2标准差 0.5 是个比较宽的设定。实际项目里我经常把 gamma 的先验收窄到 1/6 附近因为恢复率有独立文献支撑收窄它能显著改善 beta 的收敛性。R-hat 诊断需要跑 4 条不同初值的链比较链间方差与链内方差R-hat 1.05 才叫收敛只看一条链是黑匣子玄学必须多链。提示如果 accept 率低于 10%先减小 step 到 0.01 重跑如果后验呈明显的香蕉形说明 beta 和 gamma 不可辨识回第 5 章找解决办法不要硬调链长。4. 新冠疫情数据清洗与状态重建从每日新增到模型输入4.1 三类典型数据脏来源口径变化、报告延迟、周末效应任何 SDE 参数估计模型的输出质量都不会超过输入数据。新冠数据的脏我归纳成三类每一类都会直接污染参数。第一类是口径变化。疫情前半程多个地区修改过确诊定义有的把疑似合并进来有的新增“临床诊断”类目单日新增会突然跳出一根尖峰。模型看到这跟尖峰会把它解释成传播率突然升高于是 beta 的估计值被拉高后验区间却变窄这是最迷惑人的一种翻车。第二类是报告延迟。从感染到发病平均 5—7 天从发病到确诊上报又是 1—3 天。因此第 t 天报告的数字对应的其实是 t−d 天的感染事件。如果不做延迟对齐每日新增序列会比真实感染曲线整体右移导致早期增长率 r 被严重低估R0 也连带被低估。第三类是周末效应。相当多地区的检测量在周末下降、周一反弹每周七天呈现规律的锯齿。直接用原始日增数据做拟合似然函数会把这些周期波动当成过程噪声的一部分sigma_noise 的估计被顶高参数区间整体膨胀。4.2 数据预处理流程移动平均、延迟校正与累计病例我常用的清洗流程是先把累计病例做差分得到每日新增去掉负数再用中心化 7 日移动平均抹平周末效应然后按平均报告延迟把曲线左移最后截取从连续 3 天为正开始的早期窗口。代码和注释如下。import pandas as pd def preprocess_covid(df, date_col, case_col, delay_days6, min_cases5, tail_days30): 把原始疫情日增量清洗成模型可用的序列。 参数说明 ---- df : 含日期和累计/新增病例列的 DataFrame delay_days : 报告延迟天 min_cases : 连续超过这个值的日期才视为疫情开始 tail_days : 峰值后保留天数用于截断长尾 df[date_col] pd.to_datetime(df[date_col]) df df.sort_values(date_col).reset_index(dropTrue) # 1) 累计病例差分负值清零人为回填造成 df[daily] df[case_col].diff().clip(lower0) # 2) 中心化 7 日移动平均消周末效应 df[smoothed] df[daily].rolling(7, centerTrue, min_periods1).mean() # 3) 报告延迟左移t 天报告视为 t-delay_days 天感染 df[onset] df[smoothed].shift(-delay_days) # 4) 从连续 3 天 min_cases 的位置开始截取 valid df[onset] min_cases start df.loc[valid].index[0] if valid.any() else 0 # 5) 峰值后截断避免后期干预造成 beta 突变干扰早期估计 peak_idx df[onset].idxmax() end min(peak_idx tail_days, len(df) - 1) clean df.loc[start:end].reset_index(dropTrue) return clean逻辑说明第 1 步的 clip(lower0) 专门处理累计数回填造成的负日增这些负值如果不清理会让价格优化器误以为传染率突然变小。第 3 步的 shift(-delay_days) 是近似解法实际更好的做法是用反卷积把报告延迟分布解出来但日常数据量不够时均值移位已经是成本最低的对齐方案。第 5 步截断是容易被忽略的关键疫情中后期各种干预会让 beta 不再恒定而我们第 3 章估计的模型是恒定 beta硬把后期数据塞进去会让估计值变成一个说不清的中间值。加工完之后还需要一个状态重建步骤来确定 S0 和 I0。I0 我一般用新感染曲线的尾部平滑均值除以 γρ 近似I ≈ 日增感染 / (γρ)因为稳态下每日恢复人数等于每日新发。这个近似在前 30 天误差不超过 20%作为 MCMC 初值完全够用。4.3 参数初始化与边界约束让优化器不跑飞SDE 模拟的随机性使得优化器容易在参数空间里绕圈体面的初始化能省一半时间。我的初始化顺序是先算早期累计病例的对数增长率。取疫情开始后 5—15 天的累计病例对数线性回归得到斜率 r它就是 SIR 模型里的 r β − γ。然后固定 γ 为文献值 1/6 天β 的初值就是 r 1/6。from scipy import stats def init_from_growth_rate(cum_cases, gamma1/6): 用早期累计病例对数斜率估计 rbeta-gamma反推 beta 初值。 x np.arange(len(cum_cases)) y np.log(np.maximum(cum_cases, 1.0)) slope, _, _, _, _ stats.linregress(x, y) r slope beta_init r gamma return max(beta_init, 1e-3)边界约束我在代码里用 box 截断实现beta 限制在 (0.0, 3.0)gamma 在 (0.0, 1.0)。这个范围不是拍脑袋给的新冠的基本再生数大多落在 2—4 之间换算成 beta R0·γ取 γ1/6 时 R0 上限 4 就是 beta 上限 0.67我再留出 4 倍余量覆盖超传播事件gamma 低于 1/10 意味着感染期超过 10 天与绝大多数新冠文献矛盾。注意不要没有边界约束地让优化器自由游走。SDE 的随机性会造成似然面局部毛刺无约束搜索经常冲进 beta5、gamma0.02 这种离谱区域等你去解释结果时已经没法回头了。5. 参数估计避坑指南五个让模型翻车的真实场景5.1 MCMC 链不收敛参数冗余与不可辨识性现象四条链都跑了五万步R-hat 还是大于 1.1trace 图呈“毛毛虫”beta 和 gamma 的联合后验形成香蕉形分布。原因在只有日增感染观测、没有易感人群实时数据的情况下β 和 γ 是强负相关的。同样一条疫情曲线可以是“高 β 高 γ”也可以是“低 β 低 γ”两者在感染高峰期几乎无法区分。这是模型结构决定的不可辨识性不是链不够长。解决固定 γ只估计 β。恢复率 γ 的倒数就是感染周期新冠的感染周期有大量独立文献支撑5—10 天取 1/6 并把它当作已知常数后验立刻变得可收敛。另一个办法是改参数化先估计 R0 和感染周期再换算成 β 和 γ在 R0 与周期平面上采样相关性比 β-γ 平面弱得多。我实际项目里会选择后者因为 R0 的后验分布才是决策者要的产物少一次换算就少一次不确定性传播。5.2 R0 估计值剧烈漂移观测模型选错了现象用相同的数据、相同的 SDE 模型把观测误差从高斯改为负二项后R0 的后验均值从 2.9 掉到 2.3区间也窄了三分之一。原因每日新增病例是计数数据。疫情早期均值只有几十时高斯分布允许负值且把方差估计得一团糟后期均值上千时高斯又低估了极端日的方差。高斯误差让优化器以为 500 例的观测噪声和 50 例的一样大于是低值日被过度加权R0 估计被拉偏。解决统一改用负二项观测模型。负二项的方差是 μ μ²/φ允许计数数据的方差随均值走。φ 是离群参数先固定成 5—10再作为未知参数放进 MCMC。如果 φ 和过程噪声 σ_β 同时放开又会出现新的不可辨识我的经验是先固定 σ_β 为 0.1只估计 φ。5.3 预测区间过窄丢失了参数不确定性现象参数估计完成后用后验均值做一次 SDE 仿真画出的 95% 区间案例数在峰值的 70% 左右就逃出了预测带。原因我只做了单点仿真没有把 MCMC 采出的参数后验样本传播到模型输出里。过程噪声只能解释路径随机性参数不确定性才是区间的主要来源。后验均值处的仿真就像“用平均风速预测台风路径”每一条路径都会偏。解决做后验预测检验。从 MCMC 样本里随机抽 200 组参数每组参数独立仿真一条疫情曲线再把 200 条曲线的逐日分位数汇总成预测带。这个方法在第 6 章有代码。记住一个口径参数后验不确定性贡献的区间宽度通常是纯过程噪声的 3—5 倍不看这层阔度等于没做贝叶斯推断。5.4 Euler-Maruyama 步长敏感噪声强度与稳定性边界现象把内部子步 dt 从 0.1 改成 0.01beta 的后验均值不变但 sigma_noise 的估计涨了 40%把 dt 改成 0.5仿真里 S 出现负值直接被截断峰值提前两天。原因Euler-Maruyama 的强收敛阶只有 0.5。扩散项越大步长误差越大伪似然会把数值误差误判为过程噪声。sigma_noise 恰好是是那个吸走数值误差的参数。解决固定 dt≤0.05 天并把“缩短 dt 后参数变化是否超过 2%”列入模型验收标准。如果换 dt 参数变化大说明扩散项太强或者函数太刚考虑改用 Milstein 格式在 Euler-Maruyama 的更新式后面加一项 0.5·g·g′·((ΔW)²−Δt)对乘性噪声能显著提升稳定性。另外S 为负时 max(Sds, 0) 的截断虽然保住物理意义但它会人为注入正漂移所以需要在报告里说明有多少比例的时间步触发了截断超过 1% 就说明步长或噪声设置有问题。5.5 跨国数据不可比口径统一比参数更优先现象用同一个模型分别拟合三个国家的数据得到的 R0 分别是 2.1、3.6、4.9排序与当地防疫强度完全矛盾。原因三个国家报告率的差异远大于传播参数的差异。一个国家检测集中在重症另一个国家大量检测轻症同样的真实疫情到了每日新增序列上相差一个数量级。统计模型会把报告率 ρ 和 β 卷在一起β̂ 其实估的是 β·ρ 的混合体。解决给模型加报告率参数 ρ并用独立的血清学阳性率调查去锚定它。没有血清学数据时明确在论文里报告“估计的是有效报告病例层面的 β”不要直接称基本再生数。跨国比较时我只比较“校正报告率之后的相对变化”而不是绝对数值。另一个折中办法是拿同一个国家同一段时间的不同地区做相对比较地区间的报告偏误可以部分抵消。6. 验证估计结果的三个硬指标覆盖带、后验预测检验与敏感性分析6.1 后验预测检验模拟能不能重新长出新冠数据的形状参数估计完不算完要验证“模型是否真的像生成这组数据的机制”。后验预测检验PPC的做法是从后验里抽参数每套参数仿真一条完整疫情路径再加上观测噪声得到一组伪数据把这组伪数据的某个统计量与真实数据比较。统计量选三种日新增峰值、峰值到达时间、前后 14 天累计感染数。def posterior_predictive_check(samples, y_obs, S0, I0, population, n_sims200, rho0.3, phi5.0): rng np.random.default_rng(7) idx rng.choice(len(samples), sizen_sims, replaceFalse) sim_peak [] sim_peak_time [] sim_attack [] for i in idx: beta, gamma samples[i] y_sim seir_sde_daily_incidence(beta, gamma, y_obs, S0, I0, population) mu rho * y_sim y_obs_sim rng.negative_binomial(phi, phi / (phi mu)) sim_peak.append(y_obs_sim.max()) sim_peak_time.append(np.argmax(y_obs_sim)) sim_attack.append(y_obs_sim.sum()) return (sim_peak, sim_peak_time, sim_attack) obs_peak, obs_time, obs_attack (y_obs.max(), np.argmax(y_obs), y_obs.sum()) sim_peak, sim_time, sim_attack posterior_predictive_check(samples, y_obs, S0, I0, population) print(observed peak, obs_peak, vs, np.quantile(sim_peak, [0.05, 0.95])) print(observed peak time, obs_time, vs, np.quantile(sim_time, [0.05, 0.95]))逻辑说明负二项随机数生成用的是 rng.negative_binomial(phi, phi/(phimu))把过程噪声已有的路径波动和观测噪声叠加在一起。若真实峰值落在 90% 预测带外说明模型结构或观测假设有问题此时优先回头检查第 3 章的观测模型和第 4 章的延迟对齐。95% 覆盖率是另一个硬指标做 200 次仿真统计真实数据逐日落在仿真区间内的时间点比例目标在 90%—98% 之间。低于 90% 是区间过窄高于 98% 是模型把噪声放得太大都不健康。6.2 敏感性分析R0 不是唯一的决策参数传染病模型向决策者汇报时R0 只是第一行。真正影响封控时间、医疗容量准备的是峰值到达时间和累计感染负担。因此我用部分秩相关系数PRCC做敏感性分析把每个参数在其合理范围内拉丁超立方采样对每个采样点跑一次 SDE 仿真计算模型输出峰值、峰值时间、累计感染与每个参数的秩相关。输出变量最敏感参数次敏感参数工程含义峰值高度βρ减少接触率是控峰第一杠杆峰值时间γβ恢复越快高峰越早但越低14 天累计感染ρβ检测能力直接影响累计病例口径灭绝概率I0σ_β早期随机性决定是否暴发这张表看起来简单但它能立刻暴露哪些参数值得花资源去精确测量如果峰值时间对 γ 最敏感那么投资测感染周期比扩样本量更划算。R0 对决策来说是一个平均数峰值时间则是资源调配的时间表后者对管理者的价值大得多。6.3 一个落地习惯跑完估计先做三件事我现在每次完成一组疫情数据的参数估计不管多急都强制自己走三件套第一打印四条链的 R-hat 和有效样本量低于门槛的直接回炉第二从后验抽 200 组参数做 PPC把覆盖率数字贴到报告第一页第三把参数后验分布、SDE 仿真带、真实数据三条曲线叠在一张图里导出 PNG连同数据版本号和清洗脚本哈希一起存档。这个习惯救过我两次一次是发现数据源更新后旧参数不可复现一次是审稿人问“为什么你的区间这么窄”我直接甩出覆盖率记录。参数估计不是一次跑完就结束的事。数据清洗版本、延迟天数假设、噪声模型选择每改任何一处后验都会变。把验证步骤固化成流程比任何一次漂亮的拟合都值钱。希望帮到你。本文还有配套的精品资源点击获取