Julia中使用Turing包求解含强制函数ODE模型的BoundsError问题

问题原因及修复方案 #

1. 插值器越界触发BoundsError #

你当前直接用interpolate(IPTCC, BSpline(Linear()))生成的插值对象仅支持在1~36的区间取值，而DifferentialEquations的自适应求解器（Tsit5）在求解时会尝试对略超出tspan(1,36)的临时t值求值，直接调用w(t)、m(t)就会触发越界错误。
修复方法：使用绑定时间轴的插值+外插策略，超出时间范围时取边界值避免报错：

# 替换原有的插值代码段
t_axis = 1:1:total_weeks
wet_forcing = LinearInterpolation(t_axis, IPTCC, extrapolation_bc=Flat())
mobil_forcing = LinearInterpolation(t_axis, mobil, extrapolation_bc=Flat())

2. ODE函数状态索引错乱 #

你在定义epidemic_wildtype时解构的状态变量和导数赋值的对应关系存在错配，会导致后续计算逻辑异常：
原代码中状态解构为S, E, I, Hᵪ, Hₙ, R, D，其中Hᵪ为非危重住院、Hₙ为危重住院，但导数赋值、康复/死亡计算逻辑存在多处索引颠倒。修复后的ODE函数如下：

function epidemic_wildtype(dy, y, p, t)
    S, E, I, H_nc, H_c, R, D = y; # 重命名避免下标混淆：H_nc=非危重，H_c=危重
    β, λ, α, γ, θ_c, θ_nc, γ_c, γ_nc, δ_c, w, m = p;
    N = 67081000;

    dy[1] = -β*w(t)*m(t)*I*S/N + λ*R;  # S
    dy[2] = β*w(t)*m(t)*I*S/N - α*E;   # E
    dy[3] = α*E - (γ + θ_c + θ_nc)*I;     # I
    dy[4] = θ_nc*I - γ_nc*H_nc;           # 非危重住院
    dy[5] = θ_c*I - (γ_c + δ_c)*H_c;       # 危重住院
    dy[6] = γ*I + γ_nc*H_nc + γ_c*H_c - λ*R;  # R
    dy[7] = δ_c*H_c;                     # D
end

拟合段的观测匹配逻辑保持现有索引即可，已经和修复后的状态顺序对应。

3. ODE求解添加容错参数 #

为了避免求解不稳定干扰Turing采样，给solve添加容错参数适配自动微分场景：

predicted = solve(problem, Tsit5(), saveat=1, verbose=false, 
                  abstol=1e-6, reltol=1e-3, maxiters=1e5)
# 求解失败时直接给当前样本极低概率，跳过无效计算
if predicted.retcode != ReturnCode.Success
    Turing.@addlogprob! -Inf
    return
end

其他优化建议 #

• 你当前设置的链数为1，不需要使用MCMCThreads()，直接用sample(model, NUTS(0.65), 10000)即可，多链采样时再启用多线程参数。
• 可以先固定参数测试ODE求解正常后，再开启Turing采样，排查效率更高。

内容的提问来源于stack exchange，提问作者Vembha