Merge pull request #416 from davidanthoff/normalinversegaussian

Add Normal-inverse Gaussian distribution

doc/source/univariate.rst

@@ -285,6 +285,7 @@ List of Distributions
     - :ref:`logistic`
     - :ref:`lognormal`
     - :ref:`normal`
+    - :ref:`normalinversegaussian`
     - :ref:`pareto`
     - :ref:`rayleigh`
     - :ref:`symtriangular`
@@ -976,6 +977,25 @@ The probability density distribution of a `Normal distribution <http://en.wikipe
     std(d)            # Get the standard deviation, i.e. sig
 
 
+.. _normalinversegaussian:
+
+Normal-inverse Gaussian distribution
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The probability density distribution of a `Normal-inverse Gaussian distribution <http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution>`_ with location :math:`\mu`, tail heaviness :math:`\alpha`, asymmetry parameter :math:`\beta` and scale :math:`\delta` is
+
+.. math::
+
+   f(x; \mu, \alpha, \beta, \delta) = \frac{\alpha\delta K_1 \left(\alpha\sqrt{\delta^2 + (x - \mu)^2}\right)}{\pi \sqrt{\delta^2 + (x - \mu)^2}} \; e^{\delta \gamma + \beta (x - \mu)}
+
+:math:`K_j` denotes a modified Bessel function of the third kind.
+
+.. code-block:: julia
+
+    NormalInverseGaussian(mu, alpha, beta, delta)   # Normal-inverse Gaussian distribution with
+                                                    # location mu, tail heaviness alpha, asymmetry
+                                                    # parameter beta and scale delta
+
 .. _pareto:
 
 Pareto Distribution

src/Distributions.jl

@@ -112,6 +112,7 @@ export
     NormalCanon,
     NormalGamma,
     NormalInverseGamma,
+    NormalInverseGaussian,
     NormalInverseWishart,
     NormalWishart,
     Pareto,

src/univariate/continuous/normalinversegaussian.jl

@@ -0,0 +1,28 @@
+immutable NormalInverseGaussian <: ContinuousUnivariateDistribution
+  μ::Float64
+  α::Float64
+  β::Float64
+  δ::Float64
+
+  function NormalInverseGaussian(μ::Real, α::Real, β::Real, δ::Real)
+    new(μ, α, β, δ)
+  end
+end
+
+@distr_support NormalInverseGaussian -Inf Inf
+
+params(d::NormalInverseGaussian) = (d.μ, d.α, d.β, d.δ)
+
+mean(d::NormalInverseGaussian) = d.μ + d.δ * d.β / sqrt(d.α^2 - d.β^2)
+var(d::NormalInverseGaussian) = d.δ * d.α^2 / sqrt(d.α^2 - d.β^2)^3
+skewness(d::NormalInverseGaussian) = 3d.β / (d.α * sqrt(d.δ * sqrt(d.α^2 - d.β^2)))
+
+function pdf(d::NormalInverseGaussian, x::Float64)
+	μ, α, β, δ = params(d)
+	α * δ * besselk(1, α*sqrt(δ^2+(x-μ)^2)) / (π*sqrt(δ^2+(x-μ)^2)) * exp(δ*sqrt(α^2-β^2) + β*(x-μ))
+end
+
+function logpdf(d::NormalInverseGaussian, x::Float64)
+  μ, α, β, δ = params(d)
+  log(α*δ) + log(besselk(1, α*sqrt(δ^2+(x-μ)^2))) - log(π*sqrt(δ^2+(x-μ)^2)) + δ*sqrt(α^2-β^2) + β*(x-μ)
+end

src/univariates.jl

@@ -370,6 +370,7 @@ const continuous_distributions = [
     "noncentralt",
     "normal",
     "normalcanon",
+    "normalinversegaussian",
     "lognormal",    # LogNormal depends on Normal
     "pareto",
     "rayleigh",

test/normalinversegaussian.jl

@@ -0,0 +1,29 @@
+using Distributions
+using Base.Test
+
+d = NormalInverseGaussian(1.7, 1.8, 1.2, 2.3)
+
+@test_approx_eq_eps params(d)[1] 1.7 0.000000001
+@test_approx_eq_eps params(d)[2] 1.8 0.000000001
+@test_approx_eq_eps params(d)[3] 1.2 0.000000001
+@test_approx_eq_eps params(d)[4] 2.3 0.000000001
+
+# The solution was computed using this R code:
+# dnig(4.2, mu=1.7, alpha=1.8, beta=1.2, delta=2.3)
+@test_approx_eq_eps pdf(d, 4.2) 0.2021958 0.0000001
+
+# The solution was computed using this R code:
+# dnig(4.8, mu=1.7, alpha=1.8, beta=1.2, delta=2.3, log=TRUE)
+@test_approx_eq_eps logpdf(d, 4.8) -1.909973 0.000001
+
+# The solution was computed using this R code:
+# mean(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
+@test_approx_eq_eps mean(d) 3.757509 0.001
+
+# The solution was computed using this R code:
+# var(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
+@test_approx_eq_eps var(d) 3.085488 0.001
+
+# The solution was computed using this R code:
+# skewness(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
+@test_approx_eq_eps skewness(d) 1.138959 0.001

test/runtests.jl

@@ -12,6 +12,7 @@ tests = [
     "dirichlet",
     "mvnormal",
     "mvtdist",
+    "normalinversegaussian",
     "kolmogorov",
     "edgeworth",
     "matrix",