Loo cv to Bayesian Estimation of Differential Equations (challenge)

General Usage
,
1

Loo cv to Bayesian Estimation of Differential Equations (challenge)

applying cross validation with paretoSmooth.jl loo(model,chain)to models written in turing.jl works great, except for estimation in ODE(Bayesian Estimation of Differential Equations) , the error obtained is

ArgumentError: invalid base 10 digit ':' in ":,171"

that is to say, the method fails because the symbol is found: in its path, does anyone have an idea how to solve it?

2

What is the code? That looks like it’s just a typo.

1 Like
3

The code is the same as the Turing tutorial

using Turing
using DifferentialEquations

# Load StatsPlots for visualizations and diagnostics.
using StatsPlots

using LinearAlgebra

# Set a seed for reproducibility.
using Random
Random.seed!(14);

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    Ξ±, Ξ², Ξ³, Ξ΄ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (Ξ± - Ξ² * y) * x # prey
    du[2] = (Ξ΄ * x - Ξ³) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)

# Plot simulation.
#plot(solve(prob, Tsit5()))


sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

# Plot simulation and noisy observations.
plot(sol; alpha=0.3)
scatter!(sol.t, odedata'; color=[1 2], label="")

@model function fitlv(data, prob)
    # Prior distributions.
    Οƒ ~ InverseGamma(2, 3)
    Ξ± ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
    Ξ² ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
    Ξ³ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
    Ξ΄ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)

    # Simulate Lotka-Volterra model. 
    p = [Ξ±, Ξ², Ξ³, Ξ΄]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1)

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], Οƒ^2 * I)
    end

    return nothing
end

model = fitlv(odedata, prob)

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

using ParetoSmooth
loo(model,chain)
4

What’s the full error?

5

If we apply loo(model,chain) with other Turing models it works perfectly, for this particular example it generates the error:

ArgumentError: invalid base 10 digit ':' in ":,171"
6

Can you please share the error message.

7

the error message is

ArgumentError: invalid base 10 digit ':' in ":,171"

Thank you

8

I’m sure it prints a lot more than that? There’s no stack trace?

9

The full stacktrace is

ERROR: ArgumentError: invalid base 10 digit ':' in ":,51"
Stacktrace:
  [1] tryparse_internal(#unused#::Type{Int64}, s::SubString{String}, startpos::Int64, endpos::Int64, base_::Int64, raise::Bool)
    @ Base ./parse.jl:137
  [2] parse(::Type{Int64}, s::SubString{String}; base::Nothing)
    @ Base ./parse.jl:241
  [3] parse
    @ ./parse.jl:240 [inlined]
  [4] ind_from_string
    @ ~/.julia/packages/ParetoSmooth/GlhMN/src/TuringHelpers.jl:31 [inlined]
  [5] lt(o::Base.Order.By{ParetoSmooth.var"#ind_from_string#115", Base.Order.ForwardOrdering}, a::String, b::String)
    @ Base.Order ./ordering.jl:119
  [6] _issorted
    @ ./sort.jl:1174 [inlined]
  [7] _sort!(v::Vector{String}, a::Base.Sort.StableCheckSorted{Base.Sort.QuickerSort{Missing, Missing, Base.Sort.InsertionSortAlg}}, o::Base.Order.By{ParetoSmooth.var"#ind_from_string#115", Base.Order.ForwardOrdering}, kw::NamedTuple{(:scratch, :lo, :hi), Tuple{Nothing, Int64, Int64}})
    @ Base.Sort ./sort.jl:1087
  [8] _sort!
    @ ./sort.jl:692 [inlined]
  [9] _sort!
    @ ./sort.jl:643 [inlined]
 [10] _sort!
    @ ./sort.jl:714 [inlined]
 [11] _sort!
    @ ./sort.jl:659 [inlined]
 [12] _sort!
    @ ./sort.jl:595 [inlined]
 [13] #sort!#28
    @ ./sort.jl:1374 [inlined]
 [14] #sort#29
    @ ./sort.jl:1400 [inlined]
 [15] kwcall(::NamedTuple{(:by,), Tuple{ParetoSmooth.var"#ind_from_string#115"}}, ::typeof(sort), v::Vector{String})
    @ Base.Sort ./sort.jl:1400
 [16] pointwise_log_likelihoods(model::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, chains::Chains{Float64, AxisArrays.AxisArray{Float64, 3, Array{Float64, 3}, Tuple{AxisArrays.Axis{:iter, StepRange{Int64, Int64}}, AxisArrays.Axis{:var, Vector{Symbol}}, AxisArrays.Axis{:chain, UnitRange{Int64}}}}, Missing, NamedTuple{(:parameters, :internals), Tuple{Vector{Symbol}, Vector{Symbol}}}, NamedTuple{(:start_time, :stop_time), Tuple{Vector{Float64}, Vector{Float64}}}})
    @ ParetoSmooth ~/.julia/packages/ParetoSmooth/GlhMN/src/TuringHelpers.jl:33
 [17] psis_loo(::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, ::Chains{Float64, AxisArrays.AxisArray{Float64, 3, Array{Float64, 3}, Tuple{AxisArrays.Axis{:iter, StepRange{Int64, Int64}}, AxisArrays.Axis{:var, Vector{Symbol}}, AxisArrays.Axis{:chain, UnitRange{Int64}}}}, Missing, NamedTuple{(:parameters, :internals), Tuple{Vector{Symbol}, Vector{Symbol}}}, NamedTuple{(:start_time, :stop_time), Tuple{Vector{Float64}, Vector{Float64}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ParetoSmooth ~/.julia/packages/ParetoSmooth/GlhMN/src/TuringHelpers.jl:56
 [18] psis_loo(::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, ::Chains{Float64, AxisArrays.AxisArray{Float64, 3, Array{Float64, 3}, Tuple{AxisArrays.Axis{:iter, StepRange{Int64, Int64}}, AxisArrays.Axis{:var, Vector{Symbol}}, AxisArrays.Axis{:chain, UnitRange{Int64}}}}, Missing, NamedTuple{(:parameters, :internals), Tuple{Vector{Symbol}, Vector{Symbol}}}, NamedTuple{(:start_time, :stop_time), Tuple{Vector{Float64}, Vector{Float64}}}})
    @ ParetoSmooth ~/.julia/packages/ParetoSmooth/GlhMN/src/TuringHelpers.jl:55
 [19] loo(::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ParetoSmooth ~/.julia/packages/ParetoSmooth/GlhMN/src/LeaveOneOut.jl:85
 [20] loo(::DynamicPPL.Model{typeof(fitlv), (:data, :prob), (), (), Tuple{Matrix{Float64}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(lotka_volterra), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, Tuple{}, DynamicPPL.DefaultContext}, ::Vararg{Any})
    @ ParetoSmooth ~/.julia/packages/ParetoSmooth/GlhMN/src/LeaveOneOut.jl:84
 [21] top-level scope
    @ REPL[21]:1

Examining the stacktrace I’m pretty sure the issue is that this API doesn’t support non-vector data: ParetoSmooth.jl/TuringHelpers.jl at b1d603ecc34282bdfc54224247e11be3784407fe Β· TuringLang/ParetoSmooth.jl Β· GitHub. @ParadaCarleton is the best person to address this.

10

In such a case @sethaxen how can you apply loo via ArviZ.jl for the Turing example?

Thank you very much for answering

11

For ArviZ.jl, you’ll need to compute the likelihoods using Turing and then package them into an array for ArviZ. There’s an example of how to do this in Quickstart Β· ArviZ.jl, but I’ve adapted it here:

julia> using ArviZ

julia> log_likelihood = let
           log_likelihood = pointwise_loglikelihoods(
               model, MCMCChains.get_sections(chain, :parameters)
           )
           names = sort(collect(keys(log_likelihood)); by=x->parse(Int, split(x, ",")[2][1:end-1]))
           log_likelihood_data = getindex.(Ref(log_likelihood), names)
           (; data=cat(log_likelihood_data...; dims=3))
       end;

julia> idata = from_mcmcchains(chain; log_likelihood)
InferenceData with groups:
  > posterior
  > log_likelihood
  > sample_stats

julia> loo(idata)
1Γ—9 DataFrame
 Row β”‚ elpd_loo  se       p_loo    n_samples  n_data_points  warning  β‹―
     β”‚ Float64   Float64  Float64  Int64      Int64          Bool     β‹―
─────┼─────────────────────────────────────────────────────────────────
   1 β”‚ -411.174  9.93684  105.107       3000            101    false  β‹―
2 Likes
12

Excellent @sethaxen , thank you very much for your help
This is an alternative solution to what was proposed.