Problem with Turing.jl

General Usage
1

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!

2

How is your drift and diffusion defined?

3

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
4

@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
5

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
6

Hi as the previously answerd
In diffusion! You use zeros(K,L) that is a float I would replace with something like

....
dW[1:K, (2K+1):(2K+L)] = zeros(eltype(u),K,L)
...