简介：这份资源围绕传染病动力学建模中的随机微分方程展开，面向流行病建模、公共卫生干预评估方向的研究人员及高校研究生，帮助读者掌握环境噪声下疫情动态的数学刻画与数据驱动参数估计方法。压缩包共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/N，dI/dt = βSI/N - γI，dR/dt = γI。这套方程假设人群充分混合、参数恒定、没有外部干扰。但真实传播过程至少有三类随机性：

第一类是人口学随机性（demographic stochasticity），来源于个体接触和康复事件的离散性。当感染人数很少时，每次接触是否导致传播本质上是一个伯努利试验，感染人数少的时候这种随机性尤其显著，可能直接决定疫情是否消亡。

第二类是环境随机性（environmental stochasticity），来源于传播率 β 本身随时间波动。气候、人群行为变化、防控措施的松紧调整，都会让 β 不是一个常数。这类随机性影响的是参数本身，而不是事件计数。

第三类是观测噪声，来源于检测能力、报告延迟、无症状感染者的漏报。这类噪声不改变真实传播过程，但会严重扭曲你看到的数据。

SDE 主要处理的是前两类。常见做法是把环境随机性建模为对 β 的扰动：β(t) = β₀ + σ·ξ(t)，其中 ξ(t) 是白噪声。代入 SIR 后得到 Ito 型 SDE：

dS = -β₀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, seed=42): """ 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[t+1] = S[t] - infection * dt - noise I[t+1] = I[t] + (infection - recovery) * dt + noise R[t+1] = R[t] + recovery * dt # 边界保护：人数不能为负 S[t+1] = max(S[t+1], 0) I[t+1] = max(I[t+1], 0) R[t+1] = max(R[t+1], 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 中，保证 S+I+R 守恒。第三，边界保护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_sims=50): """ 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, seed=i) # 计算每日新增：I[t+1] - 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, axis=0) residual = data - sim_mean return np.sum(residual**2) # 假设 data 是观测的每日新增病例 # result = minimize(loss_function, x0=[0.3, 0.1, 0.05, 100], # args=(data,), method='Nelder-Mead', # options={'maxiter': 500})

这段代码的逻辑是：对每组候选参数，模拟多条路径，取平均后计算与观测数据的残差平方和。n_sims取 50 是一个折中，太少会导致损失函数本身噪声太大，优化器无法收敛；太多则计算成本高。method='Nelder-Mead'适合这种目标函数没有解析梯度的情况。

但最小二乘有两个问题。第一，它假设观测噪声是高斯分布且方差恒定，而实际疫情数据的噪声往往是非高斯的、异方差的。第二，它只给出了点估计，没有不确定性量化。你无法知道估计出来的 β 有多可靠。

3.3 粒子滤波做似然估计：给参数一个置信区间

要解决不确定性问题，需要从贝叶斯角度出发，计算参数的后验分布。粒子滤波（Particle Filter）是处理非线性、非高斯状态空间模型的标准工具。

核心思路是：把 SDE 离散化后看作状态空间模型，用一组粒子（样本）近似状态的后验分布，然后通过观测数据不断更新权重，最终得到参数的似然函数值。

def particle_filter_loglik(params, data, n_particles=1000, seed=0): """ 用粒子滤波计算给定参数下的对数似然 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, dtype=float) particles_I = np.full(n_particles, I0, dtype=float) particles_R = np.full(n_particles, R0, dtype=float) 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, size=n_particles, p=weights) 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_iter=5000, burn_in=1000): """ 简单的随机游走 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 samples

proposal_scale控制每步的跳跃幅度。太大接受率低，太小混合慢。目标接受率在 0.2 到 0.4 之间。burn_in是预热期，丢弃前面的样本，让链收敛到后验分布。

4. 避坑指南：SDE 参数估计中最容易翻车的五个地方

4.1 噪声项写错导致总人口不守恒

现象：模拟结果中 S+I+R 随时间漂移，有时增加有时减少。

原因：噪声项没有同时作用于 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。在均值附近做一阶泰勒展开，得到均值和协方差的 ODE：

dμ/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_draws=200): """ 从后验样本中抽取参数，模拟并计算统计量 """ rng = np.random.default_rng(456) indices = rng.choice(len(samples), size=n_draws, replace=False) 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, seed=idx) 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, axis=0)) print("模拟95%区间:", np.percentile(sim_stats, [2.5, 97.5], axis=0))

如果观测统计量落在模拟的 95% 区间内，说明模型没有系统性偏差。如果落在区间外，需要检查模型设定是否有问题。

5.3 一个我反复使用的习惯

每次做参数估计之前，我会先用模拟数据跑一遍完整流程：给定已知参数，生成模拟观测数据，然后跑参数估计，看能否恢复真实参数。这一步能帮你发现代码 bug、可识别性问题和先验设置问题。如果连模拟数据都恢复不了，真实数据就更不用想了。

这个习惯帮我省了很多时间。有一次我发现 β 的估计值总是偏高，后来用模拟数据测试才发现是观测模型中报告率的先验设得太窄，导致 β 被迫补偿。修正先验后，估计就正常了。

希望帮到你。

本文还有配套的精品资源，点击获取