SIR模型参数估计：使用autograd实现可微分ODE拟合与梯度下降优化-Python教程-PHP中文网

SIR模型参数估计：使用autograd实现可微分ODE拟合与梯度下降优化

本文介绍如何利用autograd对含常微分方程（ode）的sir流行病模型进行端到端可微分参数估计，通过替换scipy.odeint为autograd封装版本、规范参数传递方式，使梯度下降能正确计算β和γ的梯度。

本文介绍如何利用autograd对含常微分方程（ode）的sir流行病模型进行端到端可微分参数估计，通过替换scipy.odeint为autograd封装版本、规范参数传递方式，使梯度下降能正确计算β和γ的梯度。

在基于ODE的动态系统建模中，直接对数值求解器（如scipy.integrate.odeint）进行自动微分存在天然障碍——标准SciPy函数未定义梯度，且其内部实现不可被autograd追踪。若强行用原版odeint，autograd会尝试对NumPy数组操作进行符号追踪，最终因类型不匹配（如将autograd.wrappers.ArrayNode传入SciPy底层C函数）而抛出 ValueError: setting an array element with a sequence 错误。

解决的关键在于使用autograd提供的可微分ODE求解器：from autograd.scipy.integrate import odeint。该版本重写了odeint的前向与反向传播逻辑，将其注册为autograd的“原语（primitive）”，从而支持对初始条件、时间点及模型参数的精确梯度计算。

此外，还需注意两个实践细节：

参数必须以元组形式传入：odeint仅接受位置参数（positional args），不支持args=(beta, gamma)这样的关键字参数写法。需先将参数数组转为元组：
```
params = np.array([0.4, 0.05])
params_tuple = tuple(params)  # → (0.4, 0.05)
sol = odeint(sir, Y0, t, params_tuple)  # ✅ 正确调用
```
损失函数需保持纯函数特性：避免使用np.random、plt等副作用操作；所有计算路径必须可微（如用np.linalg.norm(...)而非np.sum((...)**2)更稳定，但两者均可）。

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

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

# SIR ODE 定义（使用 autograd.numpy）
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])

# 损失函数：L2误差
def loss(params, Y0, t, y_obs):
    beta, gamma = params
    sol = odeint(sir, Y0, t, tuple([beta, gamma]))
    return np.linalg.norm(y_obs - sol)

# 生成模拟观测数据
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_init = np.array([0.3, 0.03])
loss_grad = grad(loss, argnum=0)

# 梯度下降优化
learning_rate = 0.1
n_iter = 100
history = []
for i in range(n_iter):
    grads = loss_grad(params_init, Y0, t, y_obs)
    params_init -= learning_rate * grads
    history.append(params_init.copy())
    if i % 20 == 0:
        print(f"Iter {i}: β={params_init[0]:.4f}, γ={params_init[1]:.4f}")

print(f"Final estimate → β≈{params_init[0]:.4f}, γ≈{params_init[1]:.4f}")
print(f"True values  → β={true_beta:.4f}, γ={true_gamma:.4f}")

⚠️ 注意事项：