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?

**Error**

┌ 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

**Code**

@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.(μₖⱼ, σ²ₖⱼ))

return RTₖⱼ

end

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?