I’m trying to do a Bayesian inference to the following ODE system:
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
I’m using the following setup:
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)
with the data:
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)]
And based on the example from the Turing website
@model function fit_ice(data::AbstractMatrix, prob)
# Prior distributions
σₑᵣ ~ InverseGamma(2, 3)
Aₘ ~ Uniform(0.01,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
And when trying to run:
model2 = fit_ice(odedata, prob)
chain2 = sample(model2, NUTS(0.45), MCMCSerial(), 50, 2; progress=true)
I get the following error:
MethodError: ^(::Irrational{:ℯ}, ::ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 12}) is ambiguous.
I’m not sure if it stems from the fact that the ODE is nonlinear and has powers, because technically it should work since the numerical solver is able to solve the system.