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

Specific Domains Optimization (Mathematical)
,
1

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.

2

loop over both?

3

I have done that:

  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

4

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

5

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

6

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.

7

ok, perfect thanks a lot for the help