Hi all,

I am trying to use DifferentialEquations.jl and Turing.jl to perform parameter estimation using MCMC.

I am currently getting the error

TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 3}

I have seen other people get similar errors (e.g. 1, 2, can’t link to others because I am new ), and I have tried to use the suggested fixes, but I can’t seem to get it to work.

The problem seems to be caused by Turing’s automatic differentiation. I have tried to use a different AD backend (Zygote), but that did not help. I tried to turn off the AD used by the DifferentialEquations solver, but that did not help either.

Does anybody have any idea what could be going wrong and how I might fix it?

Here is a simplified version of my code:

```
using Turing, DifferentialEquations, LinearAlgebra
# Known Constants
observation_size = 100
Δ = 0.05
time = [Δ * i for i in 0:observation_size]
slices_per_block = 10
total_slices = 3 * slices_per_block
mid_of_block = ceil(Int64, slices_per_block / 2)
observed_nodes = [i * slices_per_block + mid_of_block for i in 0:2]
true_ks_int = [20., 30., 25.]
T_0 = Vector{Float64}(undef, total_slices)
T_0[1:slices_per_block] .= 330.
T_0[(slices_per_block + 1):(2 * slices_per_block)] .= 270.
T_0[(2 * slices_per_block + 1):(3 * slices_per_block)] .= 310.
# Unknown Constants
true_k12 = 5.
true_k23 = 4.
true_σ = 2.
true_p = [slices_per_block, true_ks_int, true_k12, true_k23]
function LSSM(dT, T, p, t)
s, ks_int, k12, k23 = p
N_s = 3 * s
# Put conduction coefficient between adjacent slices in list
ks = Vector{Any}(undef, N_s - 1) #I do "Any" here since k12, k23 might be ForwardDiff and not Float
ks[1:(s - 1)] .= ks_int[1]
ks[s] = k12
ks[(s + 1):(2 * s - 1)] .= ks_int[2]
ks[2 * s] = k23
ks[(2 * s + 1):(3 * s - 1)] .= ks_int[3]
# Conduction + Convection
dT[1] = (ks[1] * (T[2] - T[1]))
for i ∈ 2:(N_s - 1)
dT[i] = (ks[i - 1] * (T[i - 1] - T[i]) + ks[i] * (T[i + 1] - T[i]))
end
dT[N_s] = (ks[N_s - 1] * (T[N_s - 1] - T[N_s]))
end
LSSM_dynamics = ODEProblem(LSSM, T_0, (time[1], time[end]), true_p)
true_sol = solve(LSSM_dynamics, Tsit5(); saveat = Δ)
T_obs = true_sol[observed_nodes, :]
y = Array(T_obs) + true_σ * randn(size(Array(T_obs)))
@model function fit_LSSM(data, system, syst_consts)
σ ~ Gamma(1, 1)
k12 ~ Gamma(1, 1)
k23 ~ Gamma(1, 1)
s, ks_int, Δ = syst_consts
p = [s, ks_int, k12, k23]
predicted = solve(system, Tsit5(); p = p, saveat = Δ)
for i ∈ 1:length(predicted)
data[:, i] ~ MvNormal(predicted[observed_nodes, i], σ^2 * I)
end
end
model = fit_LSSM(y, LSSM_dynamics, [slices_per_block, true_ks_int, Δ]);
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 1; progress = false)
```

The stacktrace is very big and I don’t really understand what is going on…