Hi,

I have specified the following model in Turing.jl:

The code is:

*Model*

```
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary.Blocks
using DataInterpolations
using OrdinaryDiffEq
function System(input_funcs::Dict, params::Dict; name)
@named T_a_out = TimeVaryingFunction(input_funcs[:T_a])
@named Φ_hp_out = TimeVaryingFunction(input_funcs[:Φ_hp])
@named Φ_cv_out = TimeVaryingFunction(input_funcs[:Φ_cv])
@named Φ_s_out = TimeVaryingFunction(input_funcs[:Φ_s])
vars = @variables T_i(t) T_e(t) T_h(t) T_a(t) Φ_hp(t) Φ_cv(t) Φ_s(t)
pars = @parameters lC_i=params[:lC_i] lC_e=params[:lC_e] lC_h=params[:lC_h] lR_ie=params[:lR_ie] lR_ea=params[:lR_ea] lR_ih=params[:lR_ih] lA_w=params[:lA_w] lA_e=params[:lA_e]
eqs = [
# inputs
T_a ~ T_a_out.output.u
Φ_hp ~ Φ_hp_out.output.u
Φ_cv ~ Φ_cv_out.output.u
Φ_s ~ Φ_s_out.output.u
# dynamics
exp(lC_i)*D(T_i) ~ exp(lA_w)*Φ_s + exp(-lR_ih)*(T_h - T_i) + exp(-lR_ie)*(T_e - T_i)
exp(lC_e)*D(T_e) ~ exp(lA_e)*Φ_s + exp(-lR_ea)*(T_a - T_e) + exp(-lR_ie)*(T_i - T_e)
exp(lC_h)*D(T_h) ~ Φ_cv + Φ_hp + exp(-lR_ih)*(T_i - T_h)
]
ODESystem(eqs, t, vars, pars; systems = [T_a_out, Φ_hp_out, Φ_cv_out, Φ_s_out], name)
end
# input functions
T_a_linf = LinearInterpolation(T_a_data, time)
Φ_hp_linf = LinearInterpolation(Φ_hp_data, time)
Φ_cv_linf = LinearInterpolation(Φ_cv_data, time)
Φ_s_linf = LinearInterpolation(Φ_s_data, time)
input_funcs = Dict(
:T_a => T_a_linf,
:Φ_hp => Φ_hp_linf,
:Φ_cv => Φ_cv_linf,
:Φ_s => Φ_s_linf
)
@named system = System(input_funcs, params_d)
sys = structural_simplify(system)
prob = ODEProblem(sys, [T_bias, T_bias, T_bias], (0, time[end]))
sol = solve(prob, Tsit5());
```

*Turing spec*

```
using Turing
using DifferentialEquations
@model function LGDS(x::AbstractArray, prob, params::Dict, t_int::Float64)
# prior distributions
# V ~ Wishart(params[:ν_m], params[:Λ_m])
# Vinv ~ InverseWishart(params[:ν_m], params[:Λ_m])
σ ~ InverseGamma(2, 3)
lC ~ MvNormal(params[:μ_C], params[:Σ_C])
lR ~ MvNormal(params[:μ_R], params[:Σ_R])
lA ~ MvNormal(params[:μ_A], params[:Σ_A])
p = [lC..., lR..., lA...]
# simulate LGDS
# T_ode = solve(prob, Tsit5(); p=p, saveat=t_int, save_idxs=[1,3], verbose=false)
T_ode = solve(prob, Tsit5(); p=p, saveat=t_int, verbose=false)
if length(T_ode) != length(x)
return
end
# observations
for i in 1:length(x)
# x[:, i] ~ MvNormal(T_ode[i], Vinv)
# x[:, i] ~ MvNormalMeanPrecision(T_ode[i], V)
x[:, i] ~ MvNormalMeanPrecision(T_ode[i], σ^2 * I)
end
return nothing
end
# instantiate model
input_data = Array(data[!, [:T_i, :T_e, :T_h]])
model = LGDS(input_data, prob, prior_params, Δt)
# sample independent chains
n_chains = 3
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, n_chains; progress=true)
```

The Turing.jl code has two major problems:

- It’s very slow. The above turing code took about 80 minutes to finish.
- The posteriors are very wide, as you can see below.

I am just looking for some performance / efficiency tips, because the model is not performing well and since it’s also slow it makes it very hard to iterate.

(The model above is third order, I have also specified a first order model with only two parameters and for that model I was able to find a very accurate estimate.)