Hi! As part of GSOC, I’m working on porting the posteriordb benchmarking suite to Turing. This includes a number of Bayesian differential equations.

I’m really interested in putting Julia’s “best foot forwards”, so I’m making an effort to write my models in the most efficient and idiomatic ways possible. The DifferentialEquations models in particular seem like an area where it would be easy to miss potential optimizations.

I’m asking this here instead of on Slack, as I think this discussion could be broadly useful.

So, with that in mind, I’d appreciate feedback on the following attempt to implement a Bayesian Lotka-Volterra model in Diffeqs/Turing:

```
# Define Diffeq in DifferentialEquations form
function lotka_volterra!(dz, z, p, t)
alpha, beta, gamma, delta = p
u, v = z
# Evaluate differential equations.
dz[1] = (alpha - beta * v) * u # prey
dz[2] = (- gamma + delta * u) * v # predator
return nothing
end
# Define Turing Model
@model function lotka_volterra(N, ts, y_init, y)
theta ~ arraydist([
Truncated(Normal(1, 0.5), 0, Inf),
Truncated(Normal(0.05, 0.05), 0, Inf),
Truncated(Normal(1, 0.5), 0, Inf),
Truncated(Normal(0.05, 0.05), 0, Inf)
])
sigma ~ filldist(LogNormal(-1, 1), 2)
z_init ~ filldist(LogNormal(log(10), 1), 2)
# Create trajectory for this parameter set
prob = ODEProblem(lotka_volterra!, z_init, (0, ts[end]), theta)
z = solve(prob, DP5(), saveat = ts)
# Include trajectories in chain to match Stan
z1 := z[1, :]
z2 := z[2, :]
# If the solver failed, reject (taken roughly from DiffEqBayes.jl)
if length(z[1,:]) < N || any(z .< 0)
Turing.DynamicPPL.acclogp!!(__varinfo__, -Inf)
return
end
# Initial Condition Likelihood (y_init is observed)
for i in 1:2
y_init[i] ~ LogNormal(log(z_init[i]), sigma[i])
y[:, i] ~ MvLogNormal(log.(z[i, :]), sigma[i]^2 .* I)
end
# Generated Quantities:
# Todo, it's a bit finicky because the rand() calls break
# autodiff, so not yet set on how to do this inside the model.
end
```

I am constrained by the fact that this implementation should match the Stan implementation, which motivates some of the decisions here (e.g. making `theta`

an array dist so the resulting chains look the same, naming the parameters whatever they are named in the Stan model, using DP5 to match their solver, and so on.)

Some ideas:

- Should I avoid creating a new problem inside the model? Is this actually costly? If so, what’s the best way to set this up outside of the model?
- Putting aside the odd names and lack of special characters, which I’ve chosen to match Stan, are there any changes I should make to the code to make it align better with some style conventions which I’m accidentally violating?
- Is there some obvious optimization I’m missing?

A complete reproducible gist corresponding to the Julia code above is here.