I am trying to implement a model similar to Maas et al 2011 - equation 12 and 13 (details of my implementation below).
The issue I am facing is that I am not able to use NUTS with this model. I get repeated numerical instability error and the code just keeps running for hours (It feels the code is stuck in a loop). If I use PG(50), it takes 5-6 mins for the model to sample 4000 samples in 4 chains, but the samples just sticks closer to the prior distribution. Hence the posterior predictive distribution is pretty different than observed data. I tried MH() too, but it keeps sampling the same value. Hence, I thought to use NUTS to get better samples, but I keep getting below error. Please advise what may I be doing wrong?
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)
└ @ AdvancedHMC C:\Users.….julia\packages\AdvancedHMC\51xgc\src\hamiltonian.jl:47
@model model1(k,j, RTₖⱼ=missing) = begin
aⁱ ~ filldist(Uniform(0.01,0.5),1,j)
vⁱ ~ filldist(Uniform(0.01,100),1,j)
aᵖ ~ filldist(LogNormal(0,1),k,1)
vᵖ ~ filldist(LogNormal(0,1),k,1)
tᵉʳ ~ filldist(LogNormal(0,1),k)
μₖⱼ, σ²ₖⱼ = model_algo(aⁱ, aᵖ, vⁱ, vᵖ, tᵉʳ)
RTₖⱼ ~ arraydist(Normal.(μₖⱼ, σ²ₖⱼ))
The model_algo() implements below calculations to get the μₖⱼ, σ²ₖⱼ parameters
Also, is there a way I can debug and see what data/calculation while running NUTS is causing the issue?