# Normal Inverse Gaussian distribution in Turing

**URL:** https://discourse.julialang.org/t/normal-inverse-gaussian-distribution-in-turing/67831
**Category:** Probabilistic Programming
**Created:** [September 7, 2021, 8:00pm UTC](https://discourse.julialang.org/t/normal-inverse-gaussian-distribution-in-turing/67831 "2021-09-07T20:00:26Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![robertdj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertdj/32/103_2.png) [@robertdj](https://discourse.julialang.org/u/robertdj)
#### Post date: [September 7, 2021, 8:00pm UTC](https://discourse.julialang.org/t/normal-inverse-gaussian-distribution-in-turing/67831/1 "2021-09-07T20:00:26Z")

</div>

I’m trying to use Turing to do inference for a Normal-Inverse-Gaussian (NIG) distribution, but I have two problems:

1. I can only make the inference work for certain realizations.
2. The parameter restriction feels clumsy.

My attempt is to enforce that `|beta| < alpha`:

```julia
using Turing

@model function nig1(x)
    μ ~ Normal(0, 1)
    α ~ Exponential(1)
    δ ~ Exponential(1)

    β_limit = 0.9 * abs(α)
    β ~ Uniform(-β_limit, β_limit)
    β = clamp(β, -β_limit, β_limit)
 
    x ~ filldist(NormalInverseGaussian(μ, α, β, δ), length(x))
end

```

Without `clamp` I get errors trying to compute `sqrt` of a negative number (in `NormalInverseGaussian`’s constructor).  
If “0.9” is replaced with e.g. “0.95” there are numerical problems.

Turing can return a sensible output with this function – for certain realizations! Other realizations cause problems when evaluating the Bessel function in the pdf of the NIG.

```julia
import Random
Random.seed!(12)

μ = 0; α = 0.5; β = 0.2; δ = 1;
X = Distributions.NormalInverseGaussian(μ, α, β, δ)
data = rand(X, 2_000)
model = nig1(data)
chain = Turing.sample(model, NUTS(), 3_000)

```

The `sample` function succeeds, but there are lots of rejected proposals.

To avoid the pdf of the NIG I try to use the alternative expression as a mean-variance mixture:

```julia
@model function nig2(x)
    μ ~ Normal(0, 1)
    α ~ Exponential(1)
    δ ~ Exponential(1)

    β_limit = 0.95 * α
    β ~ Uniform(-β_limit, β_limit)
    β = clamp(β, -β_limit, β_limit)
    γ = sqrt(α^2 - β^2)

    z = filldist(InverseGaussian(δ/γ, δ^2), length(x))
    for i in eachindex(x)
        x[i] ~ Normal(μ + β * z, sqrt(z[i]))
    end
end

```

(For now) I’m not interested in the latent `z` variables, which is why I sample with `=` instead of `~`. I don’t know if this is smart?  
(Replacing the loop with `x ~ MvNormal(μ + β * z, LinearAlgebra.diagm(z))` doesn’t seem to work.)

This takes a _lot longer_ to run and the trace plots are very suspicious – only a single value is sampled for each parameter.