Estimate posterior distribution of initial condition(s)

Specific Domains Probabilistic programming
,
1

Hello,

I understood how the defining of prior distributions to parameter values and estimating the posterior distributions works using Turing. But, I was wondering if there is a way to estimate the posterior distribution for the initial condition(s)! For example, check this image below from a research article:

image
image1239×591 47.5 KB

Parameters like nu, sigma, gamma, etc., have a defined prior value. But the first row mentions the initial condition I(0) was also taken as a prior and the posterior has then been estimated along with for the parameters. In my model, I’ve tried to define something like I0 ~ Uniform(10, 1500) and passed it with the vector u0 containing other initial conditions, which are constants. But it doesn’t work and gives the error InexactError: Int64(161//1000)!

Can someone tell me how to implement it using Turing?

Thanks in advance!

2

This sounds like the correct solution to me.

This is the thing to debug – we’d need to see your actual code to understand what you wrote and why it’s throwing this error.

3

Thank you for the reply! Please find the code and the resulting error including the stacktrace:

#=
Section 1: Import required packages
=#

using Turing, Distributions, DifferentialEquations, Interpolations
using MCMCChains, Plots, StatsPlots
using CSV, XLSX, DataFrames
using Random
Random.seed!(18431)

#=
Section 2: Read the data file containing observation data
=#

my_data = DataFrame(XLSX.readtable("observation_data.xlsx","Sheet1"; infer_eltypes = true)...);

total_weeks = 36;   # Total number of time points
N = 67081000;       # Population

y_time = 1:1:total_weeks;               # Timepoints (weeks)

y_S = Float64.(my_data.Susceptible);    # Susceptible
y_S = y_S[1:total_weeks];

y_D = Float64.(my_data.Deceased);       # Deceased
y_D = y_D[1:total_weeks];

y_HC = Float64.(my_data.Hosp_critical); # Critical hospitalizations
y_HC = y_HC[1:total_weeks];

y_T = Float64.(my_data.Hosp_total);     # Total hospitalizations
y_T = y_T[1:total_weeks];
y_HNC = y_T - y_HC;                     # Non-critical hospitalizations

observation_data = [y_S y_D y_HC y_HNC];

#=
Section 3: Define the model and the respective parameters
=#

function epidemic_wildtype!(dy, y, p, t)
    S, E, I, Hᵪ, Hₙ, R, D = y;
    β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ = p;
    N = 67081000;

    dy[1] = -β*I*S/N + λ*R;  # S
    dy[2] = β*I*S/N - α*E;   # E
    dy[3] = α*E - (γ + θᵪ + θₙ)*I;     # I
    dy[4] = θₙ*I - γₙ*Hᵪ;               # HNC
    dy[5] = θᵪ*I - (γᵪ + δᵪ)*Hₙ;       # HC
    dy[6] = γ*I + γₙ*Hₙ + γᵪ*Hᵪ - λ*R;  # R
    dy[7] = δᵪ*Hᵪ;                     # D
end

#=
Section 4: Define the priors and the Bayesian model
=#

Turing.setadbackend(:forwarddiff)
@model function fitting_epidemic_wildtype!(observ_data)
    # Priors of model parameters
    β ~ truncated(Normal(0.65, 0.1), 0, 2)
    λ ~ truncated(Normal(0.5, 0.1), 0, 5)
    α ~ truncated(Normal(0.25, 0.1), 0.1, 0.5)
    γ ~ truncated(Normal(0.05, 0.1), 0, 5)
    γₙ ~ Uniform(0.05, 0.1)
    γᵪ ~ Uniform(0.05, 0.1)
    θₙ ~ Uniform(0.09, 0.75)
    θᵪ ~ Uniform(0.09, 0.75)
    δᵪ ~ Uniform(0.1, 0.8)

    p = [β, λ, α, γ, θᵪ, θₙ, γᵪ, γₙ, δᵪ];

    # Priors of standard deviations
    σ₁ ~ InverseGamma(1, 1) # Susceptible
    σ₂ ~ InverseGamma(1, 1) # Deceased
    σ₃ ~ InverseGamma(2, 3) # Critically hospitalized
    σ₄ ~ InverseGamma(1, 1) # Non-critically hospitalized

    # Initial conditions
    N = 67081000;
    S0 = N;
    # I0 = 100;
    I0 ~ Uniform(10, 1500)
    y0 = [S0, 0, I0, 0, 0, 0, 0];

    # Solve the model and compare with observed data
    problem = ODEProblem(epidemic_wildtype!, y0, (1, 36), p)
    predicted = solve(problem, Tsit5(), saveat=1.0)

    for i = 1:length(predicted)
        observ_data[i,1] ~ Normal(predicted[1,i], σ₁)
        observ_data[i,2] ~ Normal(predicted[7,i], σ₂)
        observ_data[i,3] ~ Normal(predicted[5,i], σ₃)
        observ_data[i,4] ~ Normal(predicted[4,i], σ₄)
    end
end

#=
Section 5: Run the model-inference system and save the chains
=#

model = fitting_epidemic_wildtype!(observation_data);
number_of_chains = 1;
chain = sample(model, NUTS(0.65), MCMCThreads(), 10000, number_of_chains);

When the exact same code is run with a constant I(0), the code runs without any error. Only when it is replaced with the distribution Uniform(10,1500) I get the following error:

┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:291
ERROR: LoadError: TaskFailedException

    nested task error: TaskFailedException
    Stacktrace:
     [1] wait
       @ .\task.jl:322 [inlined]
     [2] threading_run(func::Function)
       @ Base.Threads .\threadingconstructs.jl:34
     [3] macro expansion
       @ .\threadingconstructs.jl:93 [inlined]
     [4] macro expansion
       @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:342 [inlined]
     [5] (::AbstractMCMC.var"#31#41"{Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})()
       @ AbstractMCMC .\task.jl:411

        nested task error: InexactError: Int64(161//1000)
        Stacktrace:
          [1] Integer
            @ .\rational.jl:110 [inlined]
          [2] convert(#unused#::Type{Int64}, x::Rational{Int64})
            @ Base .\number.jl:7
          [3] OrdinaryDiffEq.Tsit5ConstantCache(T::Type, T2::Type)
            @ OrdinaryDiffEq C:\Users\Bharadwaj\.julia\packages\OrdinaryDiffEq\Zi9Zh\src\tableaus\low_order_rk_tableaus.jl:604
          [4] alg_cache(alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, u::Vector{Float64}, rate_prototype::Vector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Int64}, uprev::Vector{Float64}, uprev2::Vector{Float64}, f::ODEFunction{true, typeof(epidemic_wildtype!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Int64, dt::Int64, reltol::Float64, p::Vector{Float64}, calck::Bool, 
#unused#::Val{true})
            @ OrdinaryDiffEq C:\Users\Bharadwaj\.julia\packages\OrdinaryDiffEq\Zi9Zh\src\caches\low_order_rk_caches.jl:350
          [5] __init(prob::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, ODEFunction{true, typeof(epidemic_wildtype!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Int64, dtmin::Nothing, dtmax::Int64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
            @ OrdinaryDiffEq C:\Users\Bharadwaj\.julia\packages\OrdinaryDiffEq\Zi9Zh\src\solve.jl:295
          [6] #__solve#471
            @ C:\Users\Bharadwaj\.julia\packages\OrdinaryDiffEq\Zi9Zh\src\solve.jl:4 [inlined]
          [7] #solve_call#42
            @ C:\Users\Bharadwaj\.julia\packages\DiffEqBase\FtYIB\src\solve.jl:61 [inlined]
          [8] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, ODEFunction{true, typeof(epidemic_wildtype!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, 
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
            @ DiffEqBase C:\Users\Bharadwaj\.julia\packages\DiffEqBase\FtYIB\src\solve.jl:87
          [9] #solve#43
            @ C:\Users\Bharadwaj\.julia\packages\DiffEqBase\FtYIB\src\solve.jl:73 [inlined]
         [10] fitting_epidemic_wildtype!(__model__::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, __context__::DynamicPPL.SamplingContext{DynamicPPL.SampleFromUniform, DynamicPPL.DefaultContext, Random._GLOBAL_RNG}, observ_data::Matrix{Float64})
            @ Main c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_non_NPI_code.jl:90
         [11] macro expansion
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\model.jl:465 [inlined]
         [12] _evaluate
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\model.jl:448 [inlined]
         [13] evaluate_threadunsafe
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\model.jl:421 [inlined]
         [14] Model
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\model.jl:389 [inlined]
         [15] Model
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\model.jl:383 [inlined]
         [16] VarInfo
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\varinfo.jl:127 [inlined]
         [17] VarInfo
            @ C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\varinfo.jl:126 [inlined]
         [18] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}; resume_from::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
            @ DynamicPPL C:\Users\Bharadwaj\.julia\packages\DynamicPPL\RcfQU\src\sampler.jl:69
         [19] macro expansion
            @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:123 [inlined]
         [20] macro expansion
            @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\logging.jl:15 [inlined]
         [21] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nadapts,), Tuple{Int64}}})
            @ AbstractMCMC C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:114
         [22] #sample#40
            @ C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\hmc.jl:133 [inlined]
         [23] macro expansion
            @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:351 [inlined]
         [24] (::AbstractMCMC.var"#927#threadsfor_fun#42"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})(onethread::Bool)
            @ AbstractMCMC .\threadingconstructs.jl:81
         [25] (::AbstractMCMC.var"#927#threadsfor_fun#42"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}}, Vector{DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}}, Vector{Random._GLOBAL_RNG}})()
            @ AbstractMCMC .\threadingconstructs.jl:48
Stacktrace:
  [1] sync_end(c::Channel{Any})
    @ Base .\task.jl:369
  [2] macro expansion
    @ .\task.jl:388 [inlined]
  [3] macro expansion
    @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:320 [inlined]
  [4] macro expansion
    @ C:\Users\Bharadwaj\.julia\packages\ProgressLogging\6KXlp\src\ProgressLogging.jl:328 [inlined]
  [5] macro expansion
    @ C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\logging.jl:8 [inlined]
  [6] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ::MCMCThreads, N::Int64, nchains::Int64; progress::Bool, progressname::String, kwargs::Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}})
    @ AbstractMCMC C:\Users\Bharadwaj\.julia\packages\AbstractMCMC\BPJCW\src\sample.jl:314
  [7] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, ensemble::MCMCThreads, N::Int64, n_chains::Int64; chain_type::Type, progress::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:217
  [8] sample
    @ C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:217 [inlined]
  [9] #sample#6
    @ C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:202 [inlined]
 [10] sample
    @ C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:202 [inlined]
 [11] #sample#5
    @ C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:189 [inlined]
 [12] sample(model::DynamicPPL.Model{typeof(fitting_epidemic_wildtype!), (:observ_data,), (), (), Tuple{Matrix{Float64}}, Tuple{}, DynamicPPL.DefaultContext}, alg::NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}, ensemble::MCMCThreads, N::Int64, n_chains::Int64)
    @ Turing.Inference C:\Users\Bharadwaj\.julia\packages\Turing\uMQmD\src\inference\Inference.jl:189
 [13] top-level scope
    @ c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_non_NPI_code.jl:106
in expression starting at c:\Users\Bharadwaj\Indian Institute of Science\COVID-19 variant study - General\codes_wildtype_fitting\discourse_non_NPI_code.jl:106
4

So it’s definitively an issue of something being an integer when it should be a float. I bet the following will fix it: change

y0 = [S0, 0, I0, 0, 0, 0, 0]

to

y0 = zeros(eltype(I0), 7)
y0[1] = S0
y0[3] = I0

And maybe the time-interval to floats too, i.e. (1., 36.), though I didn’t think this was necessary? :man_shrugging:

2 Likes
5

It worked beautifully! :grinning: Thanks a lot for the help. Also, it looks like the time span need not necessarily contain floats.

1 Like