# Numerical errors in logit normal model using Turing.jl

**URL:** <https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602>\
**Category:** Probabilistic Programming\
**Tags:** question, turing\
**Created:** [November 1, 2019, 7:01pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602 "2019-11-01T19:01:50Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 1, 2019, 7:01pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/1 "2019-11-01T19:01:51Z")

</div>

Hello everyone,

This is my first post in a Julia forum! I have really been enjoying studying Turing.jl , Gen, and Julia in general. Thank you to all developers!

I am trying to fit a simple “logit normal” model, that is a typical Binomial model where the probabilities vary between trials. This variation is modeled as a normal distribution on the logit scale. I’m including a reproducible example below. When I run this in Julia 1.2.0 i get a large number of numerical errors, followed by what appears to be a single sample from the posterior (not 2000, as I had hoped). Any help would be appreciated !

```julia
ns = fill(400, 50)

# Sigmoid function (numerically stable)
sigmoid(x::T) where {T<:Real} = x >= zero(x) ? inv(exp(-x) + one(x)) : exp(x) / (exp(x) + one(x))

# variation in probabilities
ps = sigmoid.(rand(Normal(-1,0.5), 50))

xs = rand.(Binomial.(ns, ps))

@model normal_logit(n, x) = begin
    n_obs = length(n)

    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a = Vector{Real}(undef, n_obs)
    a ~ [Normal(ā, σ_a)]
    p = sigmoid.(a)
    
    for i in eachindex(n)
        x[i] ~ Binomial(n[i], p[i])
    end
end

normal_logit_samples = sample(normal_logit(xs, ps), NUTS(100, 0.65), 2000)

```

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 1, 2019, 7:57pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/2 "2019-11-01T19:57:26Z")

</div>

Hi, great to read that.

Regarding your question. The `sigmoid` function is unfortunately not numerically stable. Therefore, we provide the `BinomialLogit` distribution as a numerically stable alternative in Turing.

See: [https://github.com/TuringLang/Turing.jl/blob/03eba139853d950119821fbcac847b43b0caaf06/src/stdlib/distributions.jl#L45](https://github.com/TuringLang/Turing.jl/blob/03eba139853d950119821fbcac847b43b0caaf06/src/stdlib/distributions.jl#L45)

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 1, 2019, 8:36pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/3 "2019-11-01T20:36:29Z")

</div>

Btw. You can probably significantly improve the performance of sampling by using reverse mode AD and making the model type stable.

See this example:

```julia
@model gauss2(x, ::Type{TV}=Vector{Float64}) where {TV} = begin
            priors = TV(undef, 2)
            priors[1] ~ InverseGamma(2,3) # s
            priors[2] ~ Normal(0, sqrt(priors[1])) # m
            for i in 1:length(x)
                x[i] ~ Normal(priors[2], sqrt(priors[1]))
            end
        end

```

In your case you would replace the Array{Real} by the TV syntax used above.

Apologies for the short responses, I’m writing on the phone.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 1, 2019, 8:51pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/4 "2019-11-01T20:51:29Z")

</div>

> [@aammd](#):
>
> a = Vector{Real}(undef, n\_obs)  
> a ~ [Normal(ā, σ\_a)]

Change the above to:

```julia
a ~ MvNormal(fill(ā, n_obs), σ_a)

```

This together with using `Turing.setadbackend(:reverse_diff)` before sampling should give you a good speedup.

---

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 8, 2019, 6:29pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/5 "2019-11-08T18:29:39Z")

</div>

Thank you everyone who answered! after including all your suggestions I still receive many thousands of numerical errors. when the sampler does “finish”, there seems to be only a single effective sample (at least, the errors are on the order of `e-17`) !

Here is one of my error messages:

```julia
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfiniteθ = true
│ isfiniter = true
│ isfiniteℓπ = false
│ isfiniteℓκ = true
└ @ AdvancedHMC ~/.julia/packages/AdvancedHMC/PQWco/src/hamiltonian.jl:36

```

And here is a part of the result (incidentally, the estimates are not at all close to the values used in simulation)

```julia
│ 301 │ ā │ 1.27113 │ 1.99893e-15 │ 4.58585e-17 │ 0.0 │ 7.63052 │ 0.999474 │
│ 302 │ σ_a │ 3.66604 │ 4.44206e-15 │ 1.01908e-16 │ 0.0 │ 7.63052 │ 0.999474 │

```

Finally here is what I hope is a completely reproducible example. Any help would be greatly appreciated!

```julia
import Pkg; Pkg.activate(".")
using Distributions
using Turing

ns = fill(400, 300)

# Sigmoid function for generating data
sigmoid(x::T) where {T<:Real} = x >= zero(x) ? inv(exp(-x) + one(x)) : exp(x) / (exp(x) + one(x))

# variation in probabilities
ps = sigmoid.(rand(Normal(-2,0.8), 300))

xs = rand.(Binomial.(ns, ps))

@model normal_logit(n, ::Type{TV}=Vector{Float64}) where {TV} = begin

    n_obs = length(n)
    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a = TV(undef, n_obs)
    a ~ MvNormal(fill(ā, n_obs), σ_a)

    for i in eachindex(n)
        x[i] ~ BinomialLogit(n[i], a[i])
    end
end

Turing.setadbackend(:reverse_diff)

normal_logit_samples = sample(normal_logit(xs, ps), NUTS(100, 0.65), 2000)

```

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 7:08pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/6 "2019-11-08T19:08:27Z")

</div>

I’m struggling to understand your model. But I think at the moment it is not correctly implemented.

I think it might be like this, but I’m not 100% sure what you aim to do.

```julia
# your model
# n: number of trials for each observation
# x: observed value
@model normal_logit(n, x) = begin

    n_obs = length(n)

    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a ~ MvNormal(fill(ā, n_obs), σ_a)

    for i in eachindex(n)
        x[i] ~ BinomialLogit(n[i], a[i])
    end
end

```

and then…

```julia
sample(normal_logit(ns, xs), NUTS(100, 0.65), 2000)

```

This way the sampling works just fine.

Does this make sense to you?

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 7:38pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/7 "2019-11-08T19:38:09Z")

</div>

Alternatively you can do:

```julia
@model normal_logit(n, x) = begin

    n_obs = length(n)

    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a ~ MvNormal(fill(ā, n_obs), σ_a)

    x ~ VecBinomialLogit(n, a)
end

```

which will likely be faster. You can get some more performance by using FillArrays.jl.

```julia
using FillArrays
@model normal_logit(n, x) = begin

    n_obs = length(n)

    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a ~ MvNormal(Fill(ā, n_obs), σ_a)

    x ~ VecBinomialLogit(n, a)
end

```

which will be more memory efficient due to less allocations.

It’s a bit inconvenient that we call a multivariate version of `BinomialLogit` a `VecBinomialLogit`. Besides the doc string is pretty bad. I’ll open a PR fixing those issues.

---

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 8, 2019, 7:40pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/8 "2019-11-08T19:40:27Z")

</div>

It does! So the problem was actually with trying to make the model type stable in the first place? Are there any guidelines for when that is a good idea? That is something I only learned about here, not in the docs.

The model is not so very exotic; it is sometimes called a “logit normal” model, where the parameter for the Binomial distribution is modelled as an average constant plus some normally distributed error on the logit scale:

![Screenshot%20from%202019-11-08%2014-32-50](https://global.discourse-cdn.com/julialang/original/3X/0/f/0fa7f95b8e8b7d6f7256bff70aedb417fd36dbbc.png)

indeed, sampling works fine (without errors) for me, but seems to take _two hours_. In Stan I was able to draw three times as many samples in less than 5 min

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 7:41pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/9 "2019-11-08T19:41:49Z")

</div>

Try the last example I posted. Takes around 10 seconds on my 10 years old Macbook. 😃

---

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 8, 2019, 7:46pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/10 "2019-11-08T19:46:40Z")

</div>

Thank you very much! I admit I tried hard to use `VegBinomialLogit` and did not succeed.  
Here I am having a very similar error to my previous one:

```julia
MethodError: no method matching VecBinomialLogit(::Array{Int64,1}, ::TrackedArray{…,SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}})

```

output of `st` is

```julia
    Status `~/Documents/projects/MREs/Turing_binomial/Project.toml`
  [31c24e10] Distributions v0.21.6
  [1a297f60] + FillArrays v0.8.0
    Status `~/Documents/projects/MREs/Turing_binomial/Manifest.toml`
  [1a297f60] + FillArrays v0.8.0

```

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 7:49pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/11 "2019-11-08T19:49:04Z")

</div>

> [@aammd](#):
>
> It does! So the problem was actually with trying to make the model type stable in the first place? Are there any guidelines for when that is a good idea? That is something I only learned about here, not in the docs.

This and that you called the sampling with the `xs` and `ps` instead of `ns` and `xs`.

Type-stability always helps to make execution of an Julia program faster. This is nothing PPL related. If the model is not type-stable the whole type inference during sampling makes the computation expensive.

We are still having a performance tips page on the todo list.  
@mohamed82008 , are you still working on it? Do you need some help?

Some general tips:  
You should try to use multivariate distributions if possible as this will reduce allocation and computation costs. FillArrays often helps if you have a situation as in your model.  
If you cannot express it using a multivariate distribution, you should use the typed syntax I proposed in one of the first posts.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 7:52pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/12 "2019-11-08T19:52:31Z")

</div>

Please do the following for now. I’ll included this into the PR. Sorry for the inconvenience!

```julia
using FillArrays

"""
    MvBinomialLogit(n::Vector{<:Int}, logpitp::Vector{<:Real})

A multivariate binomial logit distribution with `n` trials and logit values of `logitp`.
"""
struct MvBinomialLogit{T<:Real, I<:Int} <: DiscreteMultivariateDistribution
    n::Vector{I}
    logitp::Vector{T}
end

Distributions.length(d::MvBinomialLogit) = length(d.n)

function Distributions.logpdf(d::MvBinomialLogit{<:Real}, ks::Vector{<:Integer})
    return sum(Turing.logpdf_binomial_logit.(d.n, d.logitp, ks))
end

@model normal_logit(n, x) = begin

    n_obs = length(n)

    ## hyperparameters for varying effectgs
    ā ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a ~ MvNormal(Fill(ā, n_obs), σ_a)

    x ~ MvBinomialLogit(n, a)
end

```

---

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 8, 2019, 8:44pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/13 "2019-11-08T20:44:08Z")

</div>

@trappmartin thank you so much for your help so far. However I can’t replicate your success so far! When I try with default settings, I get a warning, a long pause, followed by a progress bar that predicts _24 hours_

this is the warning:

```julia
lbeta(x::Real, w::Real)` is deprecated, use `(logabsbeta(x, w))[1]` instead.
│ caller = logpdf_binomial_logit(::Int64, [etc]

```

If i set `Turing.setadbackend(:reverse_diff)`, i get

```julia
Not implemented: convert tracked Tracker.TrackedReal{Float64} to tracked Float64
error(::String) at error.jl:33
convert(::Type{Tracker.TrackedReal{Float64}}, ::Tracker.TrackedReal{Tracker.TrackedReal{Float64}}) at real.jl:39

```

Could it be the versions of the packages involved?

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 9:02pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/14 "2019-11-08T21:02:07Z")

</div>

Yes, you should use reverse mode AD. Forward mode is slow if you have many parameters to infer.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 8, 2019, 9:02pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/15 "2019-11-08T21:02:46Z")

</div>

I’ll check later. I’m not on my computer atm.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 8, 2019, 9:24pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/16 "2019-11-08T21:24:37Z")

</div>

> [@trappmartin](#):
>
> We are still having a performance tips page on the todo list.  
> @mohamed82008 , are you still working on it? Do you need some help?

Yes! Please open a PR.

---

<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:** [November 8, 2019, 9:45pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/17 "2019-11-08T21:45:59Z")

</div>

> [@aammd](#):
>
> However I can’t replicate your success so far! When I try with default settings, I get a warning, a long pause, followed by a progress bar that predicts _24 hours_

I have also experienced poor performance with reverse diff. Would you mind sharing a full working version of your model once the issues are resolved?

---

<div class="post-metadata">

**Author:** ![aammd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aammd/32/10659_2.png) [@aammd](https://discourse.julialang.org/u/aammd)\
**Post date:** [November 8, 2019, 11:41pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/18 "2019-11-08T23:41:14Z")

</div>

@Christopher_Fisher Certainly! In fact the various incarnations of this model are all here in the same github repo. Here is the [latest one using FillArrays](https://github.com/aammd/turing_MRE/blob/master/using_fillarrays.jl)

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 9, 2019, 2:11pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/19 "2019-11-09T14:11:36Z")

</div>

Hi, here is an implementation that works fine on `#master` and should work on the latest release too. I had to remove the `FillArrays` bit again, don’t know why this made trouble and worked in the first place.

@mohamed82008 I think this should work with `FillArrays` too? Not sure why this is problematic.

```julia
using Turing
using StatsFuns

ns = fill(400, 300)

# variation in probabilities
ps = logistic.(rand(Normal(-2,0.8), 300))

xs = rand.(Binomial.(ns, ps))

struct MvBinomialLogit{T<:AbstractVector{<:Real}, I<:AbstractVector{<:Int}} <: DiscreteMultivariateDistribution
    n::I
    logitp::T
end

Distributions.length(d::MvBinomialLogit) = length(d.n)

function Distributions.logpdf(d::MvBinomialLogit, ks::AbstractVector{<:Int})
    return sum(Turing.logpdf_binomial_logit.(d.n, d.logitp, ks))
end

@model normal_logit(n, x) = begin

    n_obs = length(n)

    ## hyperparameters for varying effectgs
    b ~ Normal(0, 0.5)
    σ_a ~ Truncated(Exponential(3), 0, Inf)

    # varying intercepts for each observation
    a ~ MvNormal(fill(b, n_obs), σ_a)

    x ~ MvBinomialLogit(n, a)
end

Turing.setadbackend(:reverse_diff)
normal_logit_samples = sample(normal_logit(ns, xs), NUTS(1000, 0.8), 2000);

```

Not that I had to increase the number of adaptation iterations to prevent NUTS from diverging after the adaptation. I guess `1000` is unnecessary and less would work too, I didn’t play around much with it.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [November 9, 2019, 2:12pm UTC](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602/20 "2019-11-09T14:12:02Z")

</div>

Please let me know if this works for you.

[Next page](https://discourse.julialang.org/t/numerical-errors-in-logit-normal-model-using-turing-jl/30602.md?page=2)