DomainError: Starting point has non-finite density using DynamicHMC backend

[ClaudMor]

Hello,

I'm using DynamicHMC through DiffEqBayes.jl, but I'm experiencing an error I'm not able to overcome.
This issue is similar to this and this, and as such it could be related more to DynamicHMC.jl than to DiffEqBayes.jl, so I'll be happy to move it if necessary.

So my goal is to calibrate an epidemiological model using MCMC samplers. Here is a simplified version of the model.

using DynamicHMC, Turing, DifferentialEquations, DiffEqBayes, RecursiveArrayTools, TransformVariables

# Model parameters

β = 0.01# infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead

𝒫 = vcat([β, λ_R, λ_D]...)


# regional contact matrix and regional population

## regional contact matrix
regional_all_contact_matrix = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age
N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# model fixed parameters ( transitions)
δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
δ = vcat(repeat([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))


# Initial conditions 
i₀ = 0.075 # fraction of initial infected people in every age class
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
final_time = 20
𝒯 = (1.0,final_time); 
t  = collect(1.0:20.0)

# data used for calibration
data = [0,2,4,5,5,13,17,21,26,46,59,81,111,133,154,175,209,238,283,315]


function SIRD!(du,u,p,t)  
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initialize this parameter (death probability stratified by age, taken from literature)
    



    C = regional_all_contact_matrix 

    
    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
   dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1]
    
    # Force of infection
    Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]] 
    
    # System of equations
    @. dS = -Λ*S
    @. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
    @. dR = λ_R*(1-δ)*I 
    @. dD = λ_D*δ*I
    @. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
    

end;
# the problem builds and solves correctly
problem = ODEProblem(SIRD!,ℬ, 𝒯, 𝒫  )
solution = solve(problem, saveat = 1:20);

If I try to calibrate the three parameters using dynamic_inference, I get the following error:

# define priors
variable_parameters_priors = [Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0)]

bayesian_result_dynamic = dynamichmc_inference(problem, Tsit5(), t,data', variable_parameters_priors,
                                                as( Vector, asℝ₊, length(variable_parameters_priors)); 
                                                save_idxs = [4*4+1],
                                                σ_priors = fill(InverseGamma(2, 3), size(data', 1)),
                                                mcmc_kwargs = (initialization = (q = vcat(𝒫, repeat([0.0], size(data', 1) )), ),),)

┌ Info: finding initial optimum
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:81
DomainError with [-0.38685, -1.12035, 0.612873, 1.07036]:
Starting point has non-finite density.

Stacktrace:
 [1] local_acceptance_ratio at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\stepsize.jl:180 [inlined]
 [2] warmup at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:162 [inlined]
 [3] #36 at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:446 [inlined]
 [4] BottomRF at .\reduce.jl:81 [inlined]
 [5] afoldl at .\operators.jl:526 [inlined] (repeats 2 times)
 [6] _foldl_impl at .\tuple.jl:207 [inlined]
 [7] foldl_impl at .\reduce.jl:48 [inlined]
 [8] mapfoldl_impl at .\reduce.jl:44 [inlined]
 [9] #mapfoldl#204 at .\reduce.jl:160 [inlined]
 [10] #foldl#205 at .\reduce.jl:178 [inlined]
 [11] _warmup at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:444 [inlined]
 [12] #mcmc_keep_warmup#38 at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:520 [inlined]
 [13] macro expansion at C:\Users\claud\.julia\packages\UnPack\EkESO\src\UnPack.jl:100 [inlined]
 [14] #mcmc_with_warmup#39 at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:578 [inlined]
 [15] dynamichmc_inference(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Float64,1}, ::LinearAlgebra.Adjoint{Int64,Array{Int64,1}}, ::Array{Uniform{Float64},1}, ::TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}; σ_priors::Array{InverseGamma{Float64},1}, sample_u0::Bool, rng::Random._GLOBAL_RNG, num_samples::Int64, AD_gradient_kind::Val{:ForwardDiff}, save_idxs::Array{Int64,1}, solve_kwargs::Tuple{}, mcmc_kwargs::NamedTuple{(:initialization,),Tuple{NamedTuple{(:q,),Tuple{Array{Float64,1}}}}}) at C:\Users\claud\.julia\packages\DiffEqBayes\DNrGu\src\dynamichmc_inference.jl:105
 [16] top-level scope at In[32]:1
 [17] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

I don't understand why does DynamicHMC use negative starting point, when the trasformations and the priors require that they are positive.

Curiously, if I use DynamicHMC via Turing via DiffEqBayes, it properly works:

bayesian_result = @time turing_inference(problem, Tsit5(), t, data', variable_parameters_priors; 
                    likelihood_dist_priors = [InverseGamma(2, 3)], 
                    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
                    num_samples=1000, sampler = Turing.DynamicNUTS( ),
                    syms = [Turing.@varname(theta[i]) for i in 1:length(variable_parameters_priors)],
                    sample_u0 = false, save_idxs = [4*4+1], progress = false)

┌ Info: finding initial optimum
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:81
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\OK16j\src\solve.jl:482
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\BcprQ\src\integrator_interface.jl:322
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\OK16j\src\solve.jl:482
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\BcprQ\src\integrator_interface.jl:322
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\OK16j\src\solve.jl:482
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\BcprQ\src\integrator_interface.jl:322
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\OK16j\src\solve.jl:482
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\BcprQ\src\integrator_interface.jl:322
┌ Info: found initial stepsize
│   ϵ = 0.188
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:81
┌ Info: Starting MCMC
│   total_steps = 75
│   tuning = stepsize
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.0038
│   estimated_seconds_left = 0.28
│   ϵ = 0.188
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: Starting MCMC
│   total_steps = 25
│   tuning = stepsize and LinearAlgebra.Diagonal metric
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.013
│   estimated_seconds_left = 0.32
│   ϵ = 0.0185
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: adaptation finished
│   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.10260070201710987, 0.3571388924637585, 0.08644880409307552, 0.22561617504798276]
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:107
┌ Info: Starting MCMC
│   total_steps = 50
│   tuning = stepsize and LinearAlgebra.Diagonal metric
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.17
│   estimated_seconds_left = 8.5
│   ϵ = 0.0149
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: adaptation finished
│   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.10781336743515858, 1.3797565011999835, 0.09975860313638406, 0.17923365650618753]
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:107
┌ Info: Starting MCMC
│   total_steps = 100
│   tuning = stepsize and LinearAlgebra.Diagonal metric
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.027
│   estimated_seconds_left = 2.7
│   ϵ = 0.169
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: adaptation finished
│   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09091724910691205, 1.1149705850457625, 0.07812076795706102, 0.14441581425602598]
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:107
┌ Info: Starting MCMC
│   total_steps = 200
│   tuning = stepsize and LinearAlgebra.Diagonal metric
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.0054
│   estimated_seconds_left = 1.1
│   ϵ = 0.14
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 101
│   seconds_per_step = 0.085
│   estimated_seconds_left = 8.4
│   ϵ = 0.209
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: adaptation finished
│   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09459000674764693, 0.9130219147836826, 0.05878469020872842, 0.16351282352710864]
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:107
┌ Info: Starting MCMC
│   total_steps = 400
│   tuning = stepsize and LinearAlgebra.Diagonal metric
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.035
│   estimated_seconds_left = 14.0
│   ϵ = 0.221
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 101
│   seconds_per_step = 0.081
│   estimated_seconds_left = 24.0
│   ϵ = 0.329
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 201
│   seconds_per_step = 0.06
│   estimated_seconds_left = 12.0
│   ϵ = 0.301
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 301
│   seconds_per_step = 0.057
│   estimated_seconds_left = 5.6
│   ϵ = 0.284
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: adaptation finished
│   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09846607900835572, 1.359824677788105, 0.06657162402149377, 0.19551892156767317]
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:107
┌ Info: Starting MCMC
│   total_steps = 50
│   tuning = stepsize
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.027
│   estimated_seconds_left = 1.3
│   ϵ = 0.261
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: Starting MCMC
│   total_steps = 1000
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:112
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.025
│   estimated_seconds_left = 24.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 101
│   seconds_per_step = 0.081
│   estimated_seconds_left = 73.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 201
│   seconds_per_step = 0.077
│   estimated_seconds_left = 62.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 301
│   seconds_per_step = 0.086
│   estimated_seconds_left = 60.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 401
│   seconds_per_step = 0.077
│   estimated_seconds_left = 46.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 501
│   seconds_per_step = 0.081
│   estimated_seconds_left = 41.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 601
│   seconds_per_step = 0.077
│   estimated_seconds_left = 31.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 701
│   seconds_per_step = 0.083
│   estimated_seconds_left = 25.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 801
│   seconds_per_step = 0.086
│   estimated_seconds_left = 17.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
┌ Info: MCMC progress
│   step = 901
│   seconds_per_step = 0.08
│   estimated_seconds_left = 8.0
└ @ DynamicHMC C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\reporting.jl:130
170.734104 seconds (1.21 G allocations: 68.304 GiB, 5.46% gc time)
Object of type Chains, with data of type 1000×5×1 Array{Float64,3}

Iterations        = 1:1000
Thinning interval = 1
Chains            = 1
Samples per chain = 1000
internals         = lp
parameters        = theta[1], theta[2], theta[3], σ[1]

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters     mean     std  naive_se    mcse       ess   r_hat
  ──────────  ───────  ──────  ────────  ──────  ────────  ──────
    theta[1]   0.5326  0.0213    0.0007  0.0009  379.0684  0.9990
    theta[2]   0.0103  0.0103    0.0003  0.0004  305.5189  0.9991
    theta[3]   0.0003  0.0000    0.0000  0.0000  259.6280  0.9991
        σ[1]  11.5355  2.0806    0.0658  0.0953  384.7776  0.9997

Quantiles
  parameters    2.5%    25.0%    50.0%    75.0%    97.5%
  ──────────  ──────  ───────  ───────  ───────  ───────
    theta[1]  0.4930   0.5182   0.5311   0.5473   0.5759
    theta[2]  0.0003   0.0030   0.0072   0.0138   0.0392
    theta[3]  0.0003   0.0003   0.0003   0.0003   0.0004
        σ[1]  8.1851  10.0855  11.2343  12.7214  16.4303

I'm interested in using the DynamicHMC interface since it seems to allow for parameter initialization.

Thanks in advance

[ChrisRackauckas]

I don't think I'll be getting to this any time soon, but maybe @tpapp would like to chime in if he has an idea. The initialization runs an MAP first so that could be something, but I don't see how that could get to such a bad starting point so I am pretty confused by this behavior as well.

[tpapp]

I will look into it.

[pitmonticone]

Do you have any update? I would be very interested to know.

[tpapp]

Apologies for the delay (I am working through a large backlog) and thanks for the ping.

Conceptually, the q is a point in ℝⁿ (unconstrained), after the transformation. So you take your initial parameters, tag on a sigma, and then transform back; or transform back and tag on a sigma which will be transformed with asℝ₊. I chose the latter, and would do it like this:

variable_transformation = as(Vector, asℝ₊, length(variable_parameters_priors))
q0 = vcat(inverse(variable_transformation, 𝒫), fill(0.0, size(data', 1) ))
mcmc_kwargs = (initialization = (q = q0,),
               warmup_stages = default_warmup_stages(; local_optimization = nothing))

bayesian_result_dynamic = dynamichmc_inference(problem, Tsit5(), t,data', variable_parameters_priors,
                                               variable_transformation;
                                               save_idxs = [4*4+1],
                                               σ_priors = fill(InverseGamma(2, 3), size(data', 1)),
                                               mcmc_kwargs = mcmc_kwargs)

This skips the initial optimization, which fails. I imagine that this is happening because of a singularity that is far from the typical set, but breaks MAP. But this should be fine since we have an initial point.

Let me know if I can help with anything else.

[ClaudMor]

Hello @tpapp ,

thanks for answering and sorry for replying so late.
I implemented your fix, and it properly works for both simple models like the above one and more complex models, the only problem being that it throws an error when the initialization q0 is "too good".

ERROR: Reached maximum number of iterations while bisecting interval for ϵ.

I explain myself better in the following MWE ( partially taken from the first comment) that reproduces this behaviour.
So with your fix, now this works fine:

using DynamicHMC, Turing, DifferentialEquations, DiffEqBayes, RecursiveArrayTools, TransformVariables, Optim, DiffEqParamEstim, Plots

# Model parameters

β = 0.01# infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead

P = vcat([β, λ_R, λ_D]...)


# regional contact matrix and regional population

## regional contact matrix
regional_all_contact_matrix = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age
N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# model fixed parameters ( transitions)
δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
δ = vcat(repeat([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))


# Initial conditions 
i₀ = 0.075 # fraction of initial infected people in every age class
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
final_time = 20
𝒯 = (1.0,final_time); 
t  = collect(1.0:20.0)

# data used for calibration
data = [0,2,4,5,5,13,17,21,26,46,59,81,111,133,154,175,209,238,283,315]


function SIRD!(du,u,p,t)  
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initialize this parameter (death probability stratified by age, taken from literature)
    



    C = regional_all_contact_matrix 

    
    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
   dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1]
    
    # Force of infection
    Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]] 
    
    # System of equations
    @. dS = -Λ*S
    @. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
    @. dR = λ_R*(1-δ)*I 
    @. dD = λ_D*δ*I
    @. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
    

end;
# the problem builds and solves correctly
problem = ODEProblem(SIRD!,ℬ, 𝒯, P  )
solution = solve(problem, Tsit5(),  saveat = 1:20);

# build a faster problem
fast_problem = ODEProblem(modelingtoolkitize(problem),ℬ, 𝒯, P)
fast_solution = solve(problem, Tsit5(),  saveat = 1:20);

# define priors
const variable_parameters_priors = [Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0)]

# define transoformations
variable_transformation = as(Vector, asℝ₊, length(variable_parameters_priors))

# define initial values and kwargs
q0 = vcat(inverse(variable_transformation, P), fill(0.0, size(data', 1) ))
mcmc_kwargs = (initialization = (q = q0,),
               warmup_stages = default_warmup_stages(; local_optimization = nothing))

# run DynamicHMC
bayesian_result_dynamic = dynamichmc_inference(fast_problem, Tsit5(), t,data', variable_parameters_priors,
                                               variable_transformation;
                                               save_idxs = [4*4+1],
                                               σ_priors = fill(InverseGamma(2, 3), size(data', 1)),
                                               mcmc_kwargs = mcmc_kwargs, num_samples = 100)

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\F3qfC\src\integrator_interface.jl:328
┌ Info: found initial stepsize
└   ϵ = 0.00315
┌ Info: Starting MCMC
│   total_steps = 75
└   tuning = "stepsize"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.0016
│   estimated_seconds_left = 0.12
└   ϵ = 0.00315
┌ Info: Starting MCMC
│   total_steps = 25
└   tuning = "stepsize and LinearAlgebra.Diagonal metric"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.00021
│   estimated_seconds_left = 0.0051
└   ϵ = 0.0749
┌ Info: adaptation finished
└   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.40236231155083, 0.31137659812546126, 0.3492919363452165, 0.23136714942917988]
┌ Info: Starting MCMC
│   total_steps = 50
└   tuning = "stepsize and LinearAlgebra.Diagonal metric"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.00085
│   estimated_seconds_left = 0.042
└   ϵ = 0.0533
┌ Info: adaptation finished
└   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.25375311569626813, 0.6560920328559817, 0.3752689270387867, 0.31201639993723373]
┌ Info: Starting MCMC
│   total_steps = 100
└   tuning = "stepsize and LinearAlgebra.Diagonal metric"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.002
│   estimated_seconds_left = 0.2
└   ϵ = 0.0211
┌ Info: adaptation finished
└   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.23349928137563844, 0.49562480846134777, 0.6331264780774535, 0.18102869444738545]
┌ Info: Starting MCMC
│   total_steps = 200
└   tuning = "stepsize and LinearAlgebra.Diagonal metric"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.00018
│   estimated_seconds_left = 0.035
└   ϵ = 0.0242
┌ Info: MCMC progress
│   step = 101
│   seconds_per_step = 0.0031
│   estimated_seconds_left = 0.31
└   ϵ = 0.0733
┌ Info: adaptation finished
└   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.20678108418406718, 0.5267738626279858, 0.4572416283983815, 0.1447869395318676]
┌ Info: Starting MCMC
│   total_steps = 400
└   tuning = "stepsize and LinearAlgebra.Diagonal metric"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.00014
│   estimated_seconds_left = 0.055
└   ϵ = 0.0242
┌ Info: MCMC progress
│   step = 101
│   seconds_per_step = 0.0044
│   estimated_seconds_left = 1.3
└   ϵ = 0.0228
┌ Info: MCMC progress
│   step = 201
│   seconds_per_step = 0.0021
│   estimated_seconds_left = 0.43
└   ϵ = 0.0267
┌ Info: MCMC progress
│   step = 301
│   seconds_per_step = 0.0035
│   estimated_seconds_left = 0.34
└   ϵ = 0.0366
┌ Info: adaptation finished
└   adapted_kinetic_energy = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.24400473754408097, 0.7548261634208205, 0.5741360148155081, 0.20164617282342642]
┌ Info: Starting MCMC
│   total_steps = 50
└   tuning = "stepsize"
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 9.3e-5
│   estimated_seconds_left = 0.0046
└   ϵ = 0.0309
┌ Info: Starting MCMC
└   total_steps = 100
┌ Info: MCMC progress
│   step = 1
│   seconds_per_step = 0.0011
└   estimated_seconds_left = 0.11

But if I then try to better initialize the parameters using Optim:

# obtain optim-compatible data
const optim_data = convert(Array, VectorOfArray(data))

# define cost function
cost_function =  build_loss_objective(fast_problem, Tsit5(), L2Loss(t, optim_data'), maxiters=1000000, verbose=true, save_idxs=[4*4+1])


# define bounds that mimic the priors above
const lower_bounds = [0.01, 0.01, 0.01]
const upper_bounds = [1.0, 1.0, 1.0]

# run optim
result_optim = Optim.optimize(cost_function, lower_bounds, upper_bounds, P, Optim.Fminbox(BFGS()))

# get the calibrated parameters
P_optim  = result_optim.minimizer


# check that the parameters coming from optim correspond to a kind of good fit
calibrated_solution = solve(fast_problem, Tsit5(); saveat=t, u0=ℬ,  p=P_optim)
total_predicted_deaths = [slice[end] for slice in calibrated_solution.u]
p = plot(calibrated_solution.t, total_predicted_deaths)
plot!(calibrated_solution.t, data,seriestype=:scatter)
p

# initilialize DynamicHMC with the parameters values coming from optim
q0 = vcat(inverse(variable_transformation, P_optim), fill(0.0, size(data', 1) ))
mcmc_kwargs = (initialization = (q = q0,),
               warmup_stages = default_warmup_stages(; local_optimization = nothing))

I get:

bayesian_result_dynamic = dynamichmc_inference(fast_problem, Tsit5(), t,data', variable_parameters_priors,
                                               variable_transformation;
                                               save_idxs = [4*4+1],
                                               σ_priors = fill(InverseGamma(2, 3), size(data', 1)),
                                               mcmc_kwargs = mcmc_kwargs, num_samples = 10000)

┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\ZgbOP\src\solve.jl:479
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\F3qfC\src\integrator_interface.jl:322
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq C:\Users\claud\.julia\packages\OrdinaryDiffEq\ZgbOP\src\solve.jl:479
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\F3qfC\src\integrator_interface.jl:322
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\F3qfC\src\integrator_interface.jl:348
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ DiffEqBase C:\Users\claud\.julia\packages\DiffEqBase\F3qfC\src\integrator_interface.jl:342
ERROR: Reached maximum number of iterations while bisecting interval for ϵ.
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] bisect_stepsize(::InitialStepsizeSearch, ::DynamicHMC.var"#3#4"{DynamicHMC.Hamiltonian{GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}}},DynamicHMC.PhasePoint{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},Array{Float64,1}},Float64}, ::Float64, ::Float64, ::Float64, ::Float64) at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\stepsize.jl:125
 [3] find_initial_stepsize(::InitialStepsizeSearch, ::DynamicHMC.var"#3#4"{DynamicHMC.Hamiltonian{GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}}},DynamicHMC.PhasePoint{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},Array{Float64,1}},Float64}) at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\stepsize.jl:149
 [4] warmup at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:162 [inlined]
 [5] #36 at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:446 [inlined]
 [6] BottomRF at .\reduce.jl:81 [inlined]
 [7] afoldl at .\operators.jl:526 [inlined] (repeats 2 times)
 [8] _foldl_impl at .\tuple.jl:207 [inlined]
 [9] foldl_impl(::Base.BottomRF{DynamicHMC.var"#36#37"{DynamicHMC.SamplingLogDensity{Random._GLOBAL_RNG,LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}},DynamicHMC.NUTS{Val{:generalized}},LogProgressReport{Nothing}}}}, ::NamedTuple{(:init,),Tuple{Tuple{Tuple{},DynamicHMC.WarmupState{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},Nothing}}}}, ::Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}) at .\reduce.jl:48
 [10] mapfoldl_impl(::typeof(identity), ::DynamicHMC.var"#36#37"{DynamicHMC.SamplingLogDensity{Random._GLOBAL_RNG,LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}},DynamicHMC.NUTS{Val{:generalized}},LogProgressReport{Nothing}}}, ::NamedTuple{(:init,),Tuple{Tuple{Tuple{},DynamicHMC.WarmupState{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},Nothing}}}}, ::Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}) at .\reduce.jl:44
 [11] mapfoldl(::Function, ::Function, ::Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}; kw::Base.Iterators.Pairs{Symbol,Tuple{Tuple{},DynamicHMC.WarmupState{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},Nothing}},Tuple{Symbol},NamedTuple{(:init,),Tuple{Tuple{Tuple{},DynamicHMC.WarmupState{DynamicHMC.EvaluatedLogDensity{Array{Float64,1},Float64},GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},Nothing}}}}}) at .\reduce.jl:160
 [12] #foldl#205 at .\reduce.jl:178 [inlined]
 [13] _warmup at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:444 [inlined]
 [14] mcmc_keep_warmup(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}}, ::Int64; initialization::NamedTuple{(:q,),Tuple{Array{Float64,1}}}, warmup_stages::Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}, algorithm::DynamicHMC.NUTS{Val{:generalized}}, reporter::LogProgressReport{Nothing}) at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:520
 [15] macro expansion at C:\Users\claud\.julia\packages\UnPack\EkESO\src\UnPack.jl:100 [inlined]
 [16] mcmc_with_warmup(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#28#29"{LogDensityProblems.TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:parameters, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}},DiffEqBayes.DynamicHMCPosterior{Tsit5,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},LinearAlgebra.Adjoint{Int64,Array{Int64,1}},Array{Float64,1},Array{Uniform{Float64},1},Array{InverseGamma{Float64},1},Tuple{},Array{Int64,1}}}},Float64},Float64,4},1}}}, ::Int64; initialization::NamedTuple{(:q,),Tuple{Array{Float64,1}}}, warmup_stages::Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}, algorithm::DynamicHMC.NUTS{Val{:generalized}}, reporter::LogProgressReport{Nothing}) at C:\Users\claud\.julia\packages\DynamicHMC\n1LsS\src\mcmc.jl:578
 [17] dynamichmc_inference(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Float64,1}, ::LinearAlgebra.Adjoint{Int64,Array{Int64,1}}, ::Array{Uniform{Float64},1}, ::TransformVariables.ArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}; σ_priors::Array{InverseGamma{Float64},1}, sample_u0::Bool, rng::Random._GLOBAL_RNG, num_samples::Int64, AD_gradient_kind::Val{:ForwardDiff}, save_idxs::Array{Int64,1}, solve_kwargs::Tuple{}, mcmc_kwargs::NamedTuple{(:initialization, :warmup_stages),Tuple{NamedTuple{(:q,),Tuple{Array{Float64,1}}},Tuple{Nothing,InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{LinearAlgebra.Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}}}) at C:\Users\claud\.julia\packages\DiffEqBayes\d58SB\src\dynamichmc_inference.jl:105
 [18] top-level scope at REPL[119]:1

Increasing num_samples to 1000 and 10000 does not solve the issue. It also happens with more complex models.
How may I get around this?

Thank you very much

[tpapp]

@ClaudMor: I am not an expert on this model family, but are you sure parameters near 1 actually make sense? For Optim.optimize, you actually constrain to 1, so I am wondering if you want as𝕀 instead of asℝ₊ in the transformation.

Also, a quick fix like

q0 = vcat(inverse(variable_transformation, P_optim), fill(0.0, size(data', 1) )) .- 0.001

makes it work.

Cf boundary values are problematic in the Stan manual.

[tpapp]

Incidentally, I recently removed local optimization from DynamicHMC.jl (tpapp/DynamicHMC.jl#146). It adds very little to well-behaved posteriors (the adaptation restarts multiple times anyway, so the starting point should not matter much), and often creates problems like this.