Normal Inverse Gaussian distribution in Turing

Specific Domains Probabilistic programming
1

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:

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.

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:

@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.