# Estimate posterior distribution of initial condition(s)

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:

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?

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.

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
=#

@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

Stacktrace:
[1] wait
[3] macro expansion
[4] macro expansion
[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}})()

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)
[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})
[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{}}})
[6] #__solve#471
[7] #solve_call#42
[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}}})
[9] #solve#43
[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
[12] _evaluate
[14] Model
[15] Model
[16] VarInfo
[17] VarInfo
[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}}})
[19] macro expansion
[20] macro expansion
[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}}})
[22] #sample#40
[23] macro expansion
[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}})()
Stacktrace:
[1] sync_end(c::Channel{Any})
[2] macro expansion
[3] macro expansion
[4] macro expansion
[5] macro expansion
[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}}})
[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{}}})
[8] sample
[9] #sample#6
[10] sample
[11] #sample#5
[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
``````

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?

2 Likes

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

1 Like