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

**URL:** <https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773>\
**Category:** General Usage\
**Tags:** turing, paretosmooth\
**Created:** [January 10, 2023, 9:15pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773 "2023-01-10T21:15:26Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 10, 2023, 9:15pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/1 "2023-01-10T21:15:26Z")

</div>

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](https://turing.ml/dev/tutorials/10-bayesian-differential-equations/)) , the error obtained is

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2023, 9:46pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/2 "2023-01-10T21:46:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 10, 2023, 9:52pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/3 "2023-01-10T21:52:22Z")

</div>

The code is the same as the [Turing tutorial](https://turing.ml/dev/tutorials/10-bayesian-differential-equations/)

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2023, 9:58pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/4 "2023-01-10T21:58:55Z")

</div>

What’s the full error?

---

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 10, 2023, 10:02pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/5 "2023-01-10T22:02:09Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2023, 10:17pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/6 "2023-01-10T22:17:55Z")

</div>

Can you please share the error message.

---

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 10, 2023, 10:24pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/7 "2023-01-10T22:24:29Z")

</div>

the error message is

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

```

Thank you

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2023, 10:36pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/8 "2023-01-10T22:36:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [January 10, 2023, 11:07pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/9 "2023-01-10T23:07:40Z")

</div>

The full stacktrace is

> ****
>
> ```julia
> 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](https://github.com/TuringLang/ParetoSmooth.jl/blob/b1d603ecc34282bdfc54224247e11be3784407fe/src/TuringHelpers.jl#L30-L33). @ParadaCarleton is the best person to address this.

---

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 10, 2023, 11:44pm UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/10 "2023-01-10T23:44:01Z")

</div>

In such a case @sethaxen how can you apply loo via ArviZ.jl for the [Turing example](https://turing.ml/dev/tutorials/10-bayesian-differential-equations/)?

Thank you very much for answering

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [January 11, 2023, 12:08am UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/11 "2023-01-11T00:08:16Z")

</div>

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](https://julia.arviz.org/ArviZ/stable/quickstart), but I’ve adapted it here:

```julia
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 ⋯

```

---

<div class="post-metadata">

**Author:** ![RMF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmf/32/45530_2.png) [@RMF](https://discourse.julialang.org/u/RMF)\
**Post date:** [January 11, 2023, 2:06am UTC](https://discourse.julialang.org/t/loo-cv-to-bayesian-estimation-of-differential-equations-challenge/92773/12 "2023-01-11T02:06:41Z")

</div>

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