# Problem with Turing.jl

Hello everyone,
I am trying to estimate parameters of a SDE system using Turing.jl package, but I keep running into this issue and I cannot figure out how to solve it.
This is my code:

``````using Turing
using Plots
using StatsPlots
using RecursiveArrayTools
#Prepare the data
t_end = 100.0
n_points = length(observed_data)
dt = t_end / (n_points - 1)
t_eval = 0:dt:t_end
#Set up the SDE problem with initial epsilon values
initial_ε = [0.1, 0.5, 0.4, 0.7, 0.2]

noise_dim = 2*(K+L)

noise_prototype = zeros(noise_dim, noise_dim)

prob = SDEProblem(drift!, diffusion!, u0, (0,100), initial_ε,
noise_rate_prototype=noise_prototype)

@model function fit_sde(data, prob)
#Priors for epsilon parameters
alpha  ~ Uniform(0.0, 1.0)
beta ~ filldist(Uniform(0.0, 1.0), K)  # Assuming K epsilon values
#Construct p
p = vcat(alpha, beta)

#Sample from the prior
_prob = remake(prob, p=p)

#Solve the SDE ensemble

predicted = solve(_prob,  ISSEulerHeun(), tstops=t; p=p, saveat=0.1)
#Likelihood (assuming Gaussian noise on observations)
for j in eachindex(t)
data[:, j] ~ MvNormal(predicted[j], Diagonal(fill(0.01^2, 12)))
end
end
#Create the model instance
model = fit_sde(observed_data, prob)
#Sample from the posterior distribution (using NUTS sampler)
chain = sample(model, NUTS(0.65), 1000)
``````

and this is the error I keep getting:

``````MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5})

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.Sqrt2π)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
...

Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5})
@ Base ./number.jl:7
[2] setindex!
@ ./array.jl:1024 [inlined]
[3] macro expansion
@ ./multidimensional.jl:960 [inlined]
[4] macro expansion
@ ./cartesian.jl:64 [inlined]
[5] _unsafe_setindex!(::IndexLinear, ::Matrix{Float64}, ::Diagonal{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}}}, ::UnitRange{Int64}, ::UnitRange{Int64})
@ Base ./multidimensional.jl:955
[6] _setindex!
@ ./multidimensional.jl:944 [inlined]
[7] setindex!(::Matrix{Float64}, ::Diagonal{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}}}, ::UnitRange{Int64}, ::UnitRange{Int64})
@ Base ./abstractarray.jl:1396
[8] diffusion!(dW::Matrix{Float64}, u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 5}}, t::Float64)
@ Main ./In[195]:485
``````

Does anyone know how to fix this? I would really appreciate it.
Thank you!

How is your drift and diffusion defined?

The drift and the diffusion terms are a bit complxe and defined as follows:

``````function drift!(du, u, p, t)
alpha = p[1]
beta = p[2:(K+1)]

du[1:K] = mu_S(K, ch, N, u[(K+1):(2K)], u[1:K], u[(2K+L+1):(2K+2L)], alpha, NH, epsilon, a, omega)
du[(K+1):(2K)] = mu_I(K, ch, N, u[(K+1):(2K)], u[1:K], u[(2K+L+1):(2K+2L)], alpha, beta, NH, epsilon, a, omega)
du[(2K+1):(2K+L)] = mu_S_H(u[(2K+1):(2K+L)],u[(2K+L+1):(2K+2L)], N, u[1:K],  u[(K+1):(2K)], alpha, H, L, A, epsilon, a, omega)
du[(2K+L+1):(2K+2L)] = mu_I_H(u[(2K+L+1):(2K+2L)], u[(2K+1):(2K+L)], N, u[1:K],  u[(K+1):(2K)], alpha, H, L, A, nu, epsilon, a, omega)

end

function diffusion!(dW, u, p, t)

dW[1:K, 1:K] = sigma_SS(K, ch, N, u[(K+1):(2K)], u[1:K], u[(2K+L+1):(2K+2L)], alpha, NH, epsilon, a, omega)
dW[1:K, (K+1):(2K)] = sigma_SI(K, ch, N, u[(K+1):(2K)], u[1:K], u[(2K+L+1):(2K+2L)], H, alpha, NH, epsilon, a, omega)
dW[1:K, (2K+1):(2K+L)] = zeros(K,L)
dW[1:K, (2K+L+1):(2K+2L)] = zeros(K,L)

dW[(K+1):(2K), 1:K] = sigma_II(K, u[(K+1):(2K)], beta, N, NH)
dW[(K+1):(2K), (K+1):(2K)] = zeros(K, K)
dW[(K+1):(2K), (2K+1):(2K+L)] = zeros(K,L)
dW[(K+1):(2K), (2K+L+1):(2K+2L)] = zeros(K,L)

dW[(2K+1):(2K+L), 1:K] = zeros(L,K)
dW[(2K+1):(2K+L), (K+1):(2K)] = zeros(L,K)
dW[(2K+1):(2K+L), (2K+1):(2K+L)] = sigma_SS_HH(u[(2K+1):(2K+L)],u[(2K+L+1):(2K+2L)], N, u[1:K], u[(K+1):(2K)], alpha, H, L, A, epsilon, a, omega)
dW[(2K+1):(2K+L), (2K+L+1):(2K+2L)] = sigma_SI_HH(u[(2K+1):(2K+L)],u[(2K+L+1):(2K+2L)], N, u[1:K], u[(K+1):(2K)], alpha, H, L, A, epsilon, a, omega)

dW[(2K+L+1):(2K+2L), 1:K] = zeros(L,K)
dW[(2K+L+1):(2K+2L), (K+1):(2K)] = zeros(L,K)
dW[(2K+L+1):(2K+2L), (2K+1):(2K+L)] = sigma_II_HH(u[(2K+L+1):(2K+2L)], nu, H)
dW[(2K+L+1):(2K+2L), (2K+L+1):(2K+2L)] = zeros(L, L)

end``````

@ChrisRackauckas I think the problem is in the diffusion terms, all of the terms are the sqre root of a certain function, this is. an example

``````function sigma_II_tilde(j, i, beta, N, NH)

# Compute the inside of the square root
inside_sqrt = i[j] * beta[j]

# Return the square root
return(sqrt(inside_sqrt))
end
``````

I haven’t looked at your code in detail, and as far as I can see it’s calling lots of functions which aren’t defined in what you’ve posted, but the error combined with your use of untyped `zeros` makes me think that you are using `Float64` arrays somewhere where you need to be more flexible on types, consider:

``````julia> using ForwardDiff

julia> f1(x) = x^2
f1 (generic function with 1 method)

julia> function f2(x)
out = zeros(1)
out[1] = x
return only(x)^2
end
f2 (generic function with 1 method)

julia> function f3(x)
out = zeros(typeof(x), 1)
out[1] = x
return only(x)^2
end
f3 (generic function with 1 method)

julia> ForwardDiff.derivative(f1, 5)
10

julia> ForwardDiff.derivative(f2, 5)
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f2), Int64}, Int64, 1})

julia> ForwardDiff.derivative(f3, 5)
10
``````
1 Like

``````....