Parameter estimation with ODEProblem/Turing error

Specific Domains Modelling & Simulations
1

Hi, Trying to do parameter estimation and running into an error that I am not finding from a web search:

type or paste code hereERROR: MethodError: no method matching translate!!(::Core.SSAValue, ::Tuple{Symbol, Symbol}, ::Vector{Any}, ::Bool, ::Core.CodeInfo)
Closest candidates are:
 translate!!(::Union{Core.SSAValue, Core.SlotNumber}, ::Core.NewvarNode, ::Vector{Any}, ::Bool, ::Any) at ~/.julia/packages/Libtask/9mVaO/src/tapedfunction.jl:343
 translate!!(::Union{Core.SSAValue, Core.SlotNumber}, ::GlobalRef, ::Vector{Any}, ::Bool, ::Any) at ~/.julia/packages/Libtask/9mVaO/src/tapedfunction.jl:349
 translate!!(::Union{Core.SSAValue, Core.SlotNumber}, ::Core.SlotNumber, ::Vector{Any}, ::Bool, ::Any) at ~/.julia/packages/Libtask/9mVaO/src/tapedfunction.jl:359

So after fiddling, I managed to replicate the error in a toy model using a simple logistic map that tries to estimate initial condition y0 and rate parameter r. It seems to work OK using NUTS() but when I try SMC() the above error is encountered. In a more complex model, unfortunately NUTS also fails with the same error.

Would anyone know what I am doing incorrectly?

Please and thanks.

using DifferentialEquations, LinearAlgebra
using StatsBase, Distributions

p = (1.0, 0.2  )  # r, y0
tspan = (0.0, 10.0)

function basiclogistic!( du, u, p, t )
  du[1] = p[1] * u[1] * (1.0 - u[1] )  
end

prob = ODEProblem( basiclogistic!, [p[2]],  tspan, p )
res =  solve( prob, Tsit5(), saveat=0.1 )

# sample from the above and add some noise 
keep = sample(1:size(res)[2], 20, replace = false, ordered=true)
testdata = Array(res[keep] )' + (0.05 * randn(20))
times = res.t[keep]

# Plots.scatter(times, testdata)  #looks ok...

# estimate params
using Turing

@model function basiclogist(y, prob)
  # priors
  r ~ Normal( 1.0, 0.5 )
  y0 ~ Beta(1.0, 1.0 )
  sigma ~ truncated( Cauchy(0.0, 0.5), 1e-5, 0.5 )
  p = [r, y0]
  odesol = solve(prob, Tsit5(); p=p, saveat=0.1 )
  #likelihood
  for i in 1:length(times)
    j = findall(t -> t==times[i], odesol.t)
    if length(j) > 0
      y[i] ~ Normal( odesol.u[j[1]][1] , sigma )
    end
  end
end

mlogist = basiclogist(testdata, prob)

s = sample(mlogist, NUTS(0.8), 5 ) # joy!, though y0 is a bit off
s = sample(mlogist, SMC(), 5 ) # no joy
2

A partial solution for posterity:

It seems the tutorials online seem to be needing a slight modification (see below). Type inference in Turing seems to require the ODE parameter set to be of type “T”. SMC still not happy but NUTS() is good as well as MH().


@model function basiclogist(y, prob, ::Type{T}=Float64 ) where T
  r ~ Normal( 1.0, 0.5 )
  y0 ~ Beta(1.0, 1.0 )
  sigma ~ truncated( Cauchy(0.0, 0.5), 1e-5, 0.5 )
  p = T[r, y0]
  odesol = solve(prob, Tsit5(); p=p, saveat=0.1 )
  if odesol.retcode !== :Success
    Turing.@addlogprob! -Inf
    return nothing
  end

  for i in 1:length(times)
    j = findall(t -> t==times[i], odesol.t)
    if length(j) > 0
      y[i] ~ Normal( odesol.u[j[1]][1] , sigma )
    end
  end
end
1 Like
3

This should now auto-convert on DiffEqBase v6.91+