# Parameter Estimation 10 dimensional ODE using Turing.jl for 0/1 data

**URL:** <https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115>\
**Category:** Optimization (Mathematical)\
**Tags:** ode, turing\
**Created:** [June 3, 2024, 11:42am UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115 "2024-06-03T11:42:44Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [June 3, 2024, 11:42am UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/1 "2024-06-03T11:42:44Z")

</div>

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:

```julia

# 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:

```julia
# 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2024, 11:50am UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/2 "2024-06-03T11:50:17Z")

</div>

loop over both?

---

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [June 3, 2024, 1:33pm UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/3 "2024-06-03T13:33:23Z")

</div>

I have done that:

```julia
  for i in 1:length(t_obs)
      predicted_value = predicted(t_obs[i])
      for j in 1:length(predicted_value)
          data[j, i] ~ Bernoulli(predicted_value[j]) # Ensure it's a probability
      end
    end

```

instead of the code before. But it gives me another error:  
ERROR: MethodError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer

Which I do not know where it comes from. I have run the loop with a particular values of the parameters and just printing the data[i,j] and Bernoulli(predicted\_value[j]) and in works. So I am lost

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 4, 2024, 12:48pm UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/4 "2024-06-04T12:48:36Z")

</div>

That looks like a Turing bug. I’d open an issue.

---

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [June 4, 2024, 1:01pm UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/5 "2024-06-04T13:01:10Z")

</div>

Oh. That’s bad. Do you know how to do the same with other Julia packages?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 4, 2024, 1:03pm UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/6 "2024-06-04T13:03:24Z")

</div>

Someone will know how to make Turing work for this. But if you isolate this to something that doesn’t have an ODE solver in it you may get more help. Replace the ODE solve with just a simple analytical expression and see if you can recreate the issue.

---

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [June 4, 2024, 1:45pm UTC](https://discourse.julialang.org/t/parameter-estimation-10-dimensional-ode-using-turing-jl-for-0-1-data/115115/7 "2024-06-04T13:45:01Z")

</div>

ok, perfect thanks a lot for the help