# Turing.jl - NUTS gets stuck in "The current proposal will be rejected... isfinite.((θ, r, ℓπ, ℓκ)) = (true, true, false, true)"

**URL:** https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428
**Category:** Probabilistic Programming
**Tags:** turing
**Created:** [May 21, 2022, 7:17pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428 "2022-05-21T19:17:58Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![rkmalaiya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rkmalaiya/32/31236_2.png) [@rkmalaiya](https://discourse.julialang.org/u/rkmalaiya)
#### Post date: [May 21, 2022, 7:17pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/1 "2022-05-21T19:17:58Z")

</div>

I am trying to implement a model similar to [Maas et al 2011](https://pubmed.ncbi.nlm.nih.gov/21401290/) - 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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/4/94c6a335b77e3df4d8aaddcc80541c277a891124.png)

Also, is there a way I can debug and see what data/calculation while running NUTS is causing the issue?

---

<div class="post-metadata">

### Author: ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [May 21, 2022, 7:43pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/2 "2022-05-21T19:43:26Z")

</div>

Sometimes NUTS runs into problems if the model does not fit the data well. Are you using simulated data from the model or empirical data?

Another potential issue is the range of v:

`vⁱ ~ filldist(Uniform(0.01,100),1,j)`

Would it be reasonable to use a smaller range?

One of the simplest ways to hone in on the problem is print the parameter values inside the model block.

---

<div class="post-metadata">

### Author: ![rkmalaiya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rkmalaiya/32/31236_2.png) [@rkmalaiya](https://discourse.julialang.org/u/rkmalaiya)
#### Post date: [May 22, 2022, 4:23pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/3 "2022-05-22T16:23:06Z")

</div>

Thanks @Christopher_Fisher for the response. Here are couple of things I tried

Earlier I was using the data from Maas et al 2011. But now I tested the model using the sample generated by the model. Also, I changed my vⁱ prior to reduce the range. I also reduced the number of samples to 1 and printed the variables to see the error.

> sample\_data = Qdiffusiontest(5,2,missing)()  
> sample(Qdiffusiontest(5,2,sample\_data), NUTS(1,0.65), 1);

> ```
> vⁱ ~ filldist(Uniform(0.2,1),1,j)
> 
> ```

The error shows that σ²ₖⱼ is becoming zero. I think that must mean (based on equation 2 in my implementation) that my var(RT) becomes 0 which could be due to aᵖ or vⁱ being sampled as 0.

I printed the logs and found that **aᵖ is being sampled as 0s and NaNs.**  
My aᵖ prior is:

> aᵖ ~ filldist(LogNormal(0,1),k,1)

To trying fixing the 0 issue, I changed the σ²ₖⱼ formula to log (1.001 + …) and issue seems to be fixed for 0 case. But I am not sure why NaNs are being sampled for aᵖ and how to fix that? Please advise.

**Error logs for aᵖ only:**

> aᵖ\*\*\*\*\*\*\*\*\*\*\*\*\*\*  
> ForwardDiff.Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple  
> …  
> AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}}(0.386165161519008, **0.0,0.0,0.0,0.0** ,0.386165161519008,0.0,0.0,0.0,0.0,0.0);  
> …  
> Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple{Dy  
> …  
> (0.04060390596442086, **0.0,0.0,0.0,0.0,0.0** ,0.04060390596442086,0.0,0.0,0.0,0.0); Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple{Dy  
> …  
> (0.5975391355984172, **0.0,0.0,0.0,0.0,0.0,0.0** ,0.5975391355984172,0.0,0.0,0.0); Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple{Dyn  
> …  
> (Inf,NaN,NaN,NaN,NaN,NaN,NaN,NaN,Inf,NaN,NaN); Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple{Dyn  
> …  
> ( **Inf,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN** ,Inf,NaN);;]

> ERROR: DomainError with Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:aⁱ, :vⁱ, :aᵖ, :vᵖ, :tᵉʳ), Tuple{Dyn  
> …  
> (NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN):

> Normal: the condition σ \>= zero(σ) is not satisfied.

---

<div class="post-metadata">

### Author: ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [May 22, 2022, 5:23pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/4 "2022-05-22T17:23:03Z")

</div>

Indeed, that its very strange. I wonder what happens if you reject samples of **aᵖ** which are close to zero?

The documentation shows that you can reject a sample as follows:

```julia
using Turing
using LinearAlgebra

@model function demo(x)
    m ~ MvNormal(zero(x), I)
    if dot(m, x) < 0
        Turing.@addlogprob! -Inf
        # Exit the model evaluation early
        return
    end

    x ~ MvNormal(m, I)
    return
end

```

My suspicion is that if the sampler gets stuck at the boundary you set for rejecting samples (say .0001) , there might be challenging trade-offs due to the ratio of parameters in the definitions of the expected and variance of RT.

---

<div class="post-metadata">

### Author: ![rkmalaiya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rkmalaiya/32/31236_2.png) [@rkmalaiya](https://discourse.julialang.org/u/rkmalaiya)
#### Post date: [May 24, 2022, 2:40pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/5 "2022-05-24T14:40:53Z")

</div>

Thanks @Christopher_Fisher for the suggestion. I tried it out and it helps in the sense that the code is now sampling, but it’s sampling pretty slowly. Like really slow, may be it’s rejecting a lot of samples. Also, it’s sampling for the data that I generated from the model itself (prior predictive). But for my observed data, it’s still gets stuck in “The current proposal will be rejected…” error. I made sure that my observed data don’t have any NaN / missing.

I implemented the same code in PyMC and it sampled successfully using NUTS conditioned on my observed data. The posterior stats look great. When I tried to figure out why, I realized that [PyMC transformed](https://docs.pymc.io/en/v3/pymc-examples/examples/pymc3_howto/api_quickstart.html#Automatic-transforms-of-bounded-RVs) all my variables to have (-inf, inf) bounds and hence the sampling happened successfully. If I switch off this transform feature, I start to get similar errors as Turing NUTS.

So, is there a tutorial that I can refer to see how transformation can happen in Turing. I think it’s [Bijectors.jl](https://github.com/TuringLang/Bijectors.jl), but I have no experience in how transformation works, hence I am not sure if I am thinking about this correctly. If possible, please point me to a documentation / tutorial where I can understand how to do the transformation (similar to PyMC) in Turing.jl.

Thanks for all the help!

---

<div class="post-metadata">

### Author: ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [May 24, 2022, 9:17pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/6 "2022-05-24T21:17:04Z")

</div>

@rkmalaiya, thanks for reporting back. The transformation PyMC3 uses is interesting. I would expect that sampling would speed up if the transformation makes the posterior more similar to a multivariate normal. A while back, I used Bijectors.jl outside of Turing. However, I’m not sure how to use it in Turing for your specific use case.

Another factor is the AD backend. ReverseDiff is currently the fastest once the number of parameters is about 10 or greater (which is probably true for an IRT type of model). If you use the [Multivariate Normal](https://turing.ml/dev/docs/using-turing/performancetips) as the likelihood with ReverseDiff.jl, you should get a significant speed up.

---

<div class="post-metadata">

### Author: ![neuro\_enthusiast](https://avatars.discourse-cdn.com/v4/letter/n/43a26b/32.png) [@neuro\_enthusiast](https://discourse.julialang.org/u/neuro_enthusiast)
#### Post date: [November 22, 2022, 8:20pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/7 "2022-11-22T20:20:56Z")

</div>

Would you mind providing an intuition of why NUTS would stall if the model does not fit the data well?

---

<div class="post-metadata">

### Author: ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 22, 2022, 8:37pm UTC](https://discourse.julialang.org/t/turing-jl-nuts-gets-stuck-in-the-current-proposal-will-be-rejected-isfinite-r-true-true-false-true/81428/8 "2022-11-22T20:37:22Z")

</div>

NUTS sets up Hamiltonian dynamics for a particle using the density as a potential. One thing that can happen is that the differential equations become stiff and then oscillate out of control unless the time steps are very very small.
