Loo cv to Bayesian Estimation of Differential Equations (challenge)

RMF · January 10, 2023, 9:15pm

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?

ChrisRackauckas · January 10, 2023, 9:46pm

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

RMF · January 10, 2023, 9:52pm

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)

ChrisRackauckas · January 10, 2023, 9:58pm

What’s the full error?

RMF · January 10, 2023, 10:02pm

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"

ChrisRackauckas · January 10, 2023, 10:17pm

Can you please share the error message.

RMF · January 10, 2023, 10:24pm

the error message is

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

Thank you

ChrisRackauckas · January 10, 2023, 10:36pm

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

sethaxen · January 10, 2023, 11:07pm

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.

RMF · January 10, 2023, 11:44pm

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

Thank you very much for answering

sethaxen · January 11, 2023, 12:08am

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  ⋯

RMF · January 11, 2023, 2:06am

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