使用Turing.jl执行MCMC时遭遇MethodError错误求助

参考Turing.jl的贝叶斯ODE示例构建ice系统ODE模型进行贝叶斯推理，代码如下：

ODE系统定义 #

function ice(du, u, p, t)
    εₘ, εₖ, σ, F, d, ε = u
    Aₘ, Qₘ, n, pₘ, R, T, Eₖ, Eₘ, η, dₛₛ, Fₛₛ, ε₁, ε₂, S  = p
    
    du[1] = Aₘ * (σ/F)^n * d^(-pₘ) * ℯ^(-Qₘ/(R*T))
    du[2] = (σ - Eₖ*εₖ)/η
    du[3] = Eₘ*(S - ((σ - Eₖ*εₖ)/η) - (Aₘ * (σ/F)^n * d^(-pₘ) * ℯ^(-Qₘ/(R*T))))
    du[4] = S*(Fₛₛ - F)/ε₁
    du[5] = S*(dₛₛ - d)/ε₂
    du[6] = S
end

参数、初始条件与数据生成 #

p = [0.7, 70.0, 3.5, 0.35, 8.315, 265.0, 4.0e9, 9.0e9, 1.0e11, 0.00025, 0.6, 1.0, 1.0, 1.0e-6];
u_0 = [0.0, 0.0, 0.0, 1.0, 0.001, 0.0];
tspan = (0,1000);
prob = ODEProblem(ice, u_0, tspan , p)

# 生成带噪声的观测数据
sol = solve(prob, Rosenbrock23(); saveat=0.15);
p1 = Array(sol[3,:]) + 0.8 * randn(size(Array(sol[3,:])));
p2 = Array(sol[6,:]) + 0.8 * randn(size(Array(sol[6,:])));
odedata = [reduce(hcat, p1); reduce(hcat, p2)];

贝叶斯模型定义 #

@model function fit_ice(data::AbstractMatrix, prob)
    # Prior distributions
    σₑᵣ ~ InverseGamma(2, 3)
    Aₘ ~ Uniform(0.0,10)
    Qₘ ~ Uniform(10,1e4)
    n ~ Uniform(0.5,5.0)
    pₘ ~ Uniform(0.0,1.0)
    Eₖ ~ Uniform(1e9,9e9)
    Eₘ ~ Uniform(1e9,9e9)
    η ~ Uniform(1e9,1e12)
    dₛₛ ~ Uniform(1e-5,1e-3)
    Fₛₛ ~ Uniform(0.0,1.0)
    ε₁ ~ Uniform(0.01,1.0)
    ε₂ ~ Uniform(0.01,1.0)
    S = 1e-4
    R = 8.3145
    T = 263

    # Simulate the model but save only σ and ε
    p = [Aₘ, Qₘ, n, pₘ, R, T, Eₖ, Eₘ, η, dₛₛ, Fₛₛ, ε₁, ε₂, S]
    predicted = solve(prob, Rosenbrock23(); p=p, saveat=0.15, save_idxs=[3,6])

    # Observations
     for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], σₑᵣ^2 * I)
    end

    return nothing
end
model2 = fit_ice(odedata, prob)

执行MCMC采样时：

chain2 = sample(model2, NUTS(0.45), MCMCSerial(), 50, 2; progress=true)

出现错误：

MethodError: ^(::Irrational{:ℯ}, ::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 12}) is ambiguous.

推测问题出在ODE系统的指数项，但数值积分器能正常求解该系统，寻求修复方法。

问题根源是**ℯ作为无理数类型，与ForwardDiff的Dual类型进行幂运算时存在方法歧义**，数值积分时用的是普通浮点数，因此不会触发错误。修复方法如下：

将ODE定义中所有的ℯ^x替换为exp(x)——exp()是Julia专门的指数函数，对自动微分的Dual类型有完善的支持，可消除方法歧义。

修改后的ODE系统代码：

function ice(du, u, p, t)
    εₘ, εₖ, σ, F, d, ε = u
    Aₘ, Qₘ, n, pₘ, R, T, Eₖ, Eₘ, η, dₛₛ, Fₛₛ, ε₁, ε₂, S  = p
    
    # 替换ℯ^为exp()，提取公共计算项减少重复
    exp_term = exp(-Qₘ/(R*T))
    du[1] = Aₘ * (σ/F)^n * d^(-pₘ) * exp_term
    du[2] = (σ - Eₖ*εₖ)/η
    du[3] = Eₘ*(S - ((σ - Eₖ*εₖ)/η) - (Aₘ * (σ/F)^n * d^(-pₘ) * exp_term))
    du[4] = S*(Fₛₛ - F)/ε₁
    du[5] = S*(dₛₛ - d)/ε₂
    du[6] = S
end

修改后重新运行整个流程（重新定义prob、生成数据、构建模型、采样），即可解决该方法歧义错误。

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