基于梯度下降与Autograd的SIR模型参数估计报错求助

错误原因 #

你使用的scipy.integrate.odeint是普通数值积分函数，Autograd无法自动追踪其内部计算流程，导致求导时出现ValueError: setting an array element with a sequence错误——Autograd无法将积分结果的梯度正确转换为参数的梯度。

修正步骤 #

1. 替换可微分的ODE积分函数：使用Autograd包装的autograd.scipy.integrate.odeint，它支持自动微分追踪，能正确计算积分结果对参数的梯度。
2. 统一使用Autograd工具链：确保所有数值操作基于autograd.numpy，避免混用普通Numpy/Scipy函数破坏计算图追踪。
3. 完善梯度下降流程：补充参数初始化、学习率、迭代次数等必要配置，保证梯度下降逻辑完整。

修正后的完整代码 #

import autograd
import autograd.numpy as np
import matplotlib.pyplot as plt

# 替换为Autograd兼容的odeint
from autograd.scipy.integrate import odeint
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
    # 使用Autograd兼容的odeint计算SIR解
    sol = odeint(sir, y0=Y0, t=t, args=(beta, gamma))
    # 计算均方误差（也可换回原L2范数，仅需修改err的计算方式）
    err = np.mean((y_obs - sol) ** 2)
    return err


# 生成带噪声的观测数据
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=Y0, t=t, args=(true_beta, true_gamma))
y_obs = sol_true + np.random.normal(0, 0.05, size=sol_true.shape)

# 梯度下降参数配置
beta_init, gamma_init = 0.3, 0.05  # 参数初始猜测值
params = np.array([beta_init, gamma_init])
learning_rate = 0.01
n_iterations = 1000
loss_history = []

# 构建损失函数对参数的梯度函数
loss_grad = grad(loss, argnum=0)

# 执行梯度下降
for i in range(n_iterations):
    grads = loss_grad(params, Y0, t, y_obs)
    params -= learning_rate * grads
    current_loss = loss(params, Y0, t, y_obs)
    loss_history.append(current_loss)
    if (i+1) % 100 == 0:
        print(f"Iteration {i+1}: Loss = {current_loss:.4f}, beta = {params[0]:.4f}, gamma = {params[1]:.4f}")

# 可视化结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(t, y_obs[:, 0], 'o', label='Observed S', alpha=0.5)
plt.plot(t, y_obs[:, 1], 'o', label='Observed I', alpha=0.5)
plt.plot(t, y_obs[:, 2], 'o', label='Observed R', alpha=0.5)
sol_estimated = odeint(sir, y0=Y0, t=t, args=(params[0], params[1]))
plt.plot(t, sol_estimated[:, 0], label='Estimated S')
plt.plot(t, sol_estimated[:, 1], label='Estimated I')
plt.plot(t, sol_estimated[:, 2], label='Estimated R')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Proportion')

plt.subplot(1, 2, 2)
plt.plot(loss_history)
plt.xlabel('Iteration')
plt.ylabel('Loss')
plt.title('Loss History')
plt.show()

print(f"\nTrue parameters: beta={true_beta:.4f}, gamma={true_gamma:.4f}")
print(f"Estimated parameters: beta={params[0]:.4f}, gamma={params[1]:.4f}")

额外说明 #

• Autograd的自动微分完全依赖计算图追踪，因此必须全程使用Autograd提供的工具（如autograd.numpy、autograd.scipy.integrate.odeint），不能混用普通Scipy/Numpy函数。
• 损失函数改用均方误差是为了让梯度尺度更稳定，你也可以换回原L2范数，只需将err = np.mean((y_obs - sol)**2)改为err = np.linalg.norm(y_obs - sol, 2)。
• 若梯度下降不稳定，可尝试减小学习率或增加迭代次数，初始参数的选择也会影响收敛速度。

内容的提问来源于stack exchange，提问作者Alonso Ogueda Oliva