使用 Autograd 对 SIR 模型进行梯度下降参数估计的完整教程-Python教程-PHP中文网

使用 Autograd 对 SIR 模型进行梯度下降参数估计的完整教程

本文详解如何利用 Autograd 实现 SIR 流行病模型的参数（β 和 γ）自动微分优化，重点解决 odeint 与 Autograd 兼容性问题，适用于仅掌握 NumPy 的本科生教学场景。

本文详解如何利用 autograd 实现 sir 流行病模型的参数（β 和 γ）自动微分优化，重点解决 `odeint` 与 autograd 兼容性问题，适用于仅掌握 numpy 的本科生教学场景。

在传染病建模中，SIR 模型（Susceptible–Infected–Recovered）是理解传播动力学的基础工具。其参数 β（感染率）和 γ（康复率）往往未知，需通过观测数据反演估计。手动推导损失函数对参数的梯度复杂且易错；而借助 Autograd 可实现零代码修改的自动微分——前提是正确处理其与数值积分器的兼容性。

关键挑战在于：标准 SciPy 的 odeint 不支持 Autograd 的计算图追踪。若直接使用 scipy.integrate.odeint，Autograd 会尝试对其内部 C 实现进行符号追踪，导致 ValueError: setting an array element with a sequence 等运行时错误。

✅ 正确解法是：改用 Autograd 提供的封装版 odeint，并严格遵循其调用约定：

from autograd.scipy.integrate import odeint  # ✅ 必须使用此版本

同时注意两点：

参数必须以 元组（tuple）形式传入，而非关键字参数（如 args=(beta, gamma)）；
odeint 调用必须使用位置参数顺序：odeint(func, y0, t, params_tuple)，不可混用关键字。

以下是可直接运行的完整示例（含梯度下降训练循环）：

import autograd
import autograd.numpy as np
import matplotlib.pyplot as plt
from autograd.scipy.integrate import odeint
from autograd.builtins import tuple
from autograd import grad

def sir(y, t, beta, gamma):
    S, I, R = y
    dS_dt = -beta * S * I
    dI_dt = beta * S * I - gamma * I
    dR_dt = gamma * I
    return np.array([dS_dt, dI_dt, dR_dt])

def loss(params, Y0, t, y_obs):
    beta, gamma = params
    params_tuple = tuple([beta, gamma])  # 转为元组，适配 autograd.odeint
    sol = odeint(sir, Y0, t, params_tuple)  # ✅ 位置参数调用
    return np.linalg.norm(y_obs - sol)  # L2 损失

# 生成合成观测数据（带噪声）
np.random.seed(42)
Y0 = np.array([0.95, 0.05, 0.0])
t = np.linspace(0, 30, 101)
true_beta, true_gamma = 0.5, 1/14
sol_true = odeint(sir, Y0, t, tuple([true_beta, true_gamma]))
y_obs = sol_true + np.random.normal(0, 0.05, size=sol_true.shape)

# 初始化参数 & 定义梯度
params = np.array([0.3, 0.05])  # 初始猜测
loss_grad = grad(loss, argnum=0)
learning_rate = 0.01
n_iterations = 200

# 执行梯度下降
loss_history = []
for i in range(n_iterations):
    grads = loss_grad(params, Y0, t, y_obs)
    params -= learning_rate * grads
    loss_val = loss(params, Y0, t, y_obs)
    loss_history.append(loss_val)
    if i % 50 == 0:
        print(f"Iter {i}: β={params[0]:.4f}, γ={params[1]:.4f}, loss={loss_val:.4f}")

print(f"✅ 最终估计: β ≈ {params[0]:.4f}, γ ≈ {params[1]:.4f}")
print(f"✅ 真实值: β = {true_beta}, γ = {true_gamma:.4f}")

? 重要注意事项：