Hi,

I am trying to estimate the parameters of a ODE model, the model is N equations that represent a metapopulation, where at each node you have the probability of presence at time t. I want to compare this with the observed data which are 0-1 data, 0 for absence and 1 for presence. I am using a Bernoulli distribution for the likelihood.

This is the code for the HMC:

```
# Create the ode model
p = [0.1,0.1,0.1]
prob = DifferentialEquations.ODEProblem(fun_na, u, tspan, p)
rtol=1e-10
alg= AutoVern9(Rodas4()) #For low tolerances, Rodas4()
# Load packages
using Turing
using StatsPlots
# Estimation Bayesian --------------------------------
@model function fitlv(data, prob)
# Prior distributions.
c_d_dist ~ Uniform(0, 1)
c_h_dist ~ Uniform(0, 1)
#xi_dist ~ Uniform(0, 1)
e_dist ~ Uniform(0, 1)
# Simulate Lotka-Volterra model.
p = [c_d_dist, c_h_dist, e_dist]
predicted = solve(prob, alg; p=p, reltol = rtol) # Observations.
for i in 1:length(t_obs)
data[:, i] ~ Bernoulli(predicted(t_obs[i]))
end
return nothing
end
model = fitlv(matrix_obs, prob)
# Sample posterior
iterations = 1000
ϵ = 0.05
τ = 10
chain = sample(model, Gibbs(HMC(ϵ, τ, :p), MH(:y)), iterations; progress = true)
```

This is the model:

```
# Ode function non autonomous
function fun_na(du, u, p, t)
R_M_vec = [interpolated_functions[i](t) for i in 1:N]
mat = exp.(-alp*d_1)
mat[diagind(mat)] .= 0 # Set to zero diagonal of distance
Cd = p[1].*mat*u # Natural dispersal
Ch = p[2].*(m_c.*eta)*u # Human mobility
du.= R_M_vec.*(Cd .+ Ch) .* (1 .- u) - p[3] .* u
end
```

But I think the problem is that the predicted(t_obs[I]) inside the Bernoulli it is a vector, since it is for all the nodes. But it does not understand it. How should I change it? I saw in other examples that they used the MvNormal but there is not MvBernoilli.