# Loo cv to Bayesian Estimation of Differential Equations (challenge)

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?

What is the code? That looks like itβs just a typo.

1 Like

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)
``````

Whatβs the full error?

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"
``````

Can you please share the error message.

the error message is

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

Thank you

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

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.

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

Thank you very much for answering

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

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