企业官网建设流程全解析

简介：这份资源围绕传染病动力学建模中的随机微分方程展开，面向流行病建模、公共卫生干预评估方向的研究人员及高校研究生，帮助读者掌握环境噪声下疫情传播的随机建模思路与数据驱动的参数估计方法。压缩包共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/N，dI/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, dt=0.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), size=len(state0)) states[i+1] = states[i] + drift * dt + diff * dW return states

dt 的选择很关键。dt 太大，离散化误差会累积，轨迹偏离真实 SDE 的解；dt 太小，计算量爆炸。我一般先用 dt=0.01 跑一遍，再把 dt 减半看结果是否稳定。如果两次结果差异明显，说明 dt 还不够小。Milstein 方法在扩散项依赖状态时精度更高，但需要计算扩散项的导数，实现复杂。对于 SIR 这类扩散项相对简单的模型，Euler-Maruyama 在 dt=0.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, sigma=2.0): """ 对每日报告病例做平滑处理 case_series: 每日新增报告数 sigma: 高斯平滑核宽度 """ smoothed = gaussian_filter1d(case_series.astype(float), sigma=sigma) # 保证非负 smoothed = np.maximum(smoothed, 0) return smoothed

sigma 控制平滑程度。sigma 太小，噪声没压住；sigma 太大，早期增长信号被抹平。我一般从 sigma=1 开始，逐步加到 3，看拟合残差是否还有明显结构。如果残差图呈现周期性波动，说明平滑不够；如果峰值被削平，说明平滑过度。

3.2 最小二乘估计：目标函数与优化器选择

有了平滑后的感染曲线，下一步是找一组参数 (β, γ, σ) 使模型模拟的轨迹最贴近数据。最小二乘的目标函数是模拟轨迹与观测轨迹的加权平方误差。由于 SDE 每次模拟结果不同，目标函数本身是随机的，不能直接用梯度下降。

常见做法是固定随机种子，对同一组参数多次模拟取平均，再用优化器搜索。优化器我一般用 Nelder-Mead 或差分进化，它们不需要梯度，对随机目标函数更鲁棒。

from scipy.optimize import differential_evolution def loss_function(params, observed_I, state0, t_span, n_sim=10): """ 计算参数对应的损失 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), dt=0.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))), maxiter=100, tol=1e-6)

n_sim 的选择是个权衡。n_sim 太小，损失函数噪声大，优化器容易卡在局部；n_sim 太大，计算时间成倍增加。我一般先用 n_sim=5 粗搜，找到大致范围后再用 n_sim=20 精搜。differential_evolution 的 maxiter 和 tol 控制收敛，maxiter 太小结果不稳定，太大浪费时间。tol 设 1e-6 通常够用。

3.3 贝叶斯参数估计：先验选择与后验采样

最小二乘给的是点估计，贝叶斯方法给的是后验分布。对于 SDE 模型，似然函数没有解析形式，必须用近似贝叶斯计算（ABC）或粒子滤波。ABC 的思路是：从先验采样参数，模拟轨迹，如果模拟轨迹与观测的距离小于阈值 ε，就接受这个参数样本。

def abc_rejection(observed_I, state0, t_span, prior_ranges, n_samples=10000, epsilon=0.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), dt=0.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_particles=1000, param_bounds=None): """ 粒子滤波做 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, p=weights) 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_sim=5) res = differential_evolution(loss, [(0.05, 0.5), (0.001, 0.2)], maxiter=50) profile.append(res.fun) return np.array(profile)

beta_values 的范围根据先验或粗搜结果确定，一般取 20 到 50 个点。每个点上的优化用 differential_evolution 可能太慢，可以改用 Nelder-Mead，因为只有两个参数。轮廓似然曲线如果出现多个极小，说明模型可能存在多解，需要检查数据是否覆盖了足够多的动态阶段。

我自己的习惯是：每次做 SDE 参数估计，先跑轮廓似然看参数可辨识性，再决定是否值得做完整的贝叶斯推断。如果轮廓似然平坦，贝叶斯后验也会很宽，这时候与其花时间采样，不如回头检查数据质量或简化模型。这个习惯帮我省下了不少算力，也避免了很多「后验不收敛」的玄学问题。希望帮到你。

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