I am trying to estimate the paremeters of a model using Turing, in the code below there is a toy model for simplicity. This model as output has a 10 dimensional reponse variable. Therefore for each time I have a vector as the response variable. My response is a probability of having x species at each patch, so it ranges from 0 to 1. And the observations are if in this patch there was x species so it is 1 or there is not x species. This is why I use a Bernoulli distribution for the likelihood. I have seen that for other distributions there is the Multivariate distribution such as for the normal MvNormal, but I have not seen this with the Bernoulli. So I code it myself doing two loops, as a recommedation of a julia user.

This is the code:

```
# Toy model
function fun_na(u, p, t_vec)
y = Array{Float64}(undef,length(u), length(t_vec))
for i in 1:length(t_vec)
y[:, i] = p[1].*rand(length(u))
end
return y
end
# Set parameters
t0=0.0
tf=6574.0
tspan = (t0, tf)
t_vect=1:tf
u = pop_init
simulations = fun_na(u,p,t_vect) # Example
# Load packages
using Turing
using StatsPlots
# Estimation Bayesian --------------------------------
@model function fitlv(data, t_obs)
# Prior distributions.
c_d_dist ~ Uniform(0.1e-8, 0.1)
c_h_dist ~ Uniform(0.1e-8, 0.1)
#xi_dist ~ Uniform(0, 1)
e_dist ~ Uniform(0.1e-8, 0.1)
# Simulate Lotka-Volterra model.
p = [c_d_dist, c_h_dist, e_dist]
predicted = fun_na(u,p,t_vect) # Simulations.
if length(predicted) > 0
for i in 1:length(t_obs)
predicted_value = predicted[:,i]
for j in 1:length(predicted_value)
data[j, i] ~ Bernoulli(predicted_value[j]) # Ensure it's a probability
end
end
end
# Example in Turing tutorial which does not work because it is not Multivariate
#for i in 1:length(t_obs)
# data[:, i] ~ Bernoulli(predicted[:,i])
#end
return nothing
end
# matrix_obs is a 10 x 18 matrix and t_obs are the 18 moments in time
model = fitlv(matrix_obs,t_obs)
# Sample 3 independent chains with forward-mode automatic differentiation (the default).
iterations = 1000
ϵ = 0.05
τ = 10
chain = sample(model, Gibbs(HMC(ϵ, τ, :p), MH(:y)), iterations; progress = true)
```

I get the folllowing error:

ERROR: MethodError: reducing over an empty collection is not allowed; consider supplying `init`

to the reducer

In another issue in the julia community I asked the same question with a Dynamic model much more complicated. And they have told me to use a simple toy model and ask again because it looks like a Turing bug.