I am trying to sample a difference equation using Turing but encounter errors.

Basically, it’s a SIR model that solves fine, with a quick solver function written below (I wrote my own very simple solver to understand the basics and see how the solver could be altered/extended for modeling purposes).

Then, I add random noise to the solution and use it as data for the Turing model, in order to find the parameters of the SIR model using Turing. I used this strategy successfully with the Lokta Voltera model or the tutorial (see my modification of the tutorial here), but I can’t make it work with my own simple solver.

Any help is appreciated.

## Julia version: 1.5.0.

## Error

```
chain = sample(model, NUTS(.65),1000)
ERROR: StackOverflowError:
Stacktrace:
[1] ScalMat at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:5 [inlined] (repeats 2 times)
[2] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:13 [inlined]
[3] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:14 [inlined]
[4] MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64}) at /Users/aamzallag/.julia/packages/Distributions/jFoHB/src/multivariate/mvnormal.jl:202 (repeats 79984 times)
```

## Code

```
using Turing, Distributions
# Import MCMCChain, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots
# Set a seed for reproducibility.
using Random
Random.seed!(12);
# Defining the dynamics
function sir_model(du,u,p,t)
S, I, R = u
β, γ = p
du[1] = dS = -β*S*I/N
du[2] = dI = β*S*I/N - γ*I
du[3] = dR = γ*I
end
### trying a simple discrete time solver
# this solver takes a function func_ describing
# the dynamics (a la du/dt = f(u,p,t)) and simulate
# the dynamics with a simple loop, and store the values
# in variable u
function diffsolve(func_, u0, p, n=100)
du = [0., 0., 0.]
u_t = copy(u0)
#u = [u0]
u = Array{Array{Real,1},1}() # if Real not specified,
# I get the error: in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,3}
push!(u, u0)
for t in 1:n
func_(du, u_t, p, 1)
u_t = u_t + du
push!(u, u_t)
end
#hcat(u...)
u
end
N=10e3 # total population
conf=10 # confirmed cases at the beginning
u0 = [N-conf,conf,0.0] #.* N
## Trying the solver and plotting
x=diffsolve(sir_model, u0, [0.8, 0.1], 50)
x=hcat(x...)
# add noise
#odedata = x + 0.1*maximum(x) .* randn(size(x))
odedata = x .* exp.(0.1 .* randn(size(x)))
plot(x')
scatter!(odedata')
Turing.setadbackend(:forwarddiff)
@model function fitsir(data)
σ ~ InverseGamma(1, 1)
#σ ~ truncated(Normal(0., 0.5),0,5)
β ~ truncated(Normal(0., 0.5),0,5)
γ ~ truncated(Normal(0., 0.5),0,5)
p = [β,γ]
predicted = diffsolve(sir_model, u0, p, 10)
for i = 1:length(predicted)
data[:,i] ~ MvNormal(predicted[i], σ)
end
end
model = fitsir(odedata)
chain = sample(model, NUTS(.65),1000)
# ERROR: StackOverflowError:
# Stacktrace:
# [1] ScalMat at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:5 [inlined] (repeats 2 times)
# [2] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:13 [inlined]
# [3] convert at /Users/aamzallag/.julia/packages/PDMats/8bndI/src/scalmat.jl:14 [inlined]
# [4] MvNormal(::Array{Real,1}, ::PDMats.ScalMat{Float64}) at /Users/aamzallag/.julia/packages/Distributions/jFoHB/src/multivariate/mvnormal.jl:202 (repeats 79984 times)
```