Inverse Gaussian distribution

Inverse Gaussian
Probability density function
Cumulative distribution function
Notation	${\displaystyle \textstyle \operatorname {IG} \left(\mu ,\lambda \right)}$
Parameters	${\displaystyle \textstyle \mu >0}$ ${\displaystyle \lambda >0}$
Support	${\displaystyle \textstyle x\in (0,\infty )}$
PDF	${\displaystyle \textstyle {\sqrt {\frac {\lambda }{2\pi x^{3}}}}\exp \left[-{\frac {\lambda (x-\mu )^{2}}{2\mu ^{2}x}}\right]}$
CDF	${\displaystyle \textstyle \Phi \left({\sqrt {\frac {\lambda }{x}}}\left({\frac {x}{\mu }}-1\right)\right)}$ ${\displaystyle \textstyle {}+\exp \left({\frac {2\lambda }{\mu }}\right)\Phi \left(-{\sqrt {\frac {\lambda }{x}}}\left({\frac {x}{\mu }}+1\right)\right)}$ where ${\displaystyle \Phi }$ is the standard normal (standard Gaussian) distribution c.d.f.
Mean	${\displaystyle \textstyle \operatorname {E} [X]=\mu }$ ${\displaystyle \textstyle \operatorname {E} \left[{\frac {1}{X}}\right]={\frac {1}{\mu }}+{\frac {1}{\lambda }}}$
Mode	${\displaystyle \textstyle \mu \left[\left(1+{\frac {9\mu ^{2}}{4\lambda ^{2}}}\right)^{\frac {1}{2}}-{\frac {3\mu }{2\lambda }}\right]}$
Variance	${\displaystyle \textstyle \operatorname {Var} [X]={\frac {\mu ^{3}}{\lambda }}}$ ${\displaystyle \textstyle \operatorname {Var} \left[{\frac {1}{X}}\right]={\frac {1}{\mu \lambda }}+{\frac {2}{\lambda ^{2}}}}$
Skewness	${\displaystyle \textstyle 3\left({\frac {\mu }{\lambda }}\right)^{1/2}}$
Excess kurtosis	${\displaystyle \textstyle {\frac {15\mu }{\lambda }}}$
MGF	${\displaystyle \textstyle \exp \left[{{\frac {\lambda }{\mu }}\left(1-{\sqrt {1-{\frac {2\mu ^{2}t}{\lambda }}}}\right)}\right]}$
CF	${\displaystyle \textstyle \exp \left[{{\frac {\lambda }{\mu }}\left(1-{\sqrt {1-{\frac {2\mu ^{2}\mathrm {i} t}{\lambda }}}}\right)}\right]}$

In probability theory, the inverse Gaussian distribution (also known as the Wald distribution) is a two-parameter family of continuous probability distributions with support on ⁠${\displaystyle (0,\infty )}$⁠.

Its probability density function is given by

${\displaystyle f(x;\mu ,\lambda )={\sqrt {\frac {\lambda }{2\pi x^{3}}}}\exp {\biggl (}-{\frac {\lambda (x-\mu )^{2}}{2\mu ^{2}x}}{\biggr )}}$

for ⁠${\displaystyle x>0}$⁠, where ${\displaystyle \mu >0}$ is the mean and ${\displaystyle \lambda >0}$ is a shape parameter.¹ Either ⁠${\displaystyle \mu }$⁠ or ⁠${\displaystyle \lambda }$⁠ (or more generally any combination of the form ${\displaystyle \mu ^{p}\lambda ^{1-p}}$ for any real ⁠${\displaystyle p}$⁠) can serve as a scale parameter, so a proper (i.e., unscaled) shape parameter would be any non-zero power of ⁠${\displaystyle \varphi =\lambda /\mu }$⁠: Tweedie proposed to use the ⁠${\displaystyle (\mu ,\varphi )}$⁠ and ⁠${\displaystyle (\varphi ,\lambda )}$⁠ parametrizations in addition to the standard ⁠${\displaystyle (\mu ,\lambda )}$⁠ parametrization (“Each of these forms is convenient or suggestive for some purpose.”²), and later on uses exclusively the ⁠${\displaystyle (\varphi ,\lambda )}$⁠ parametrization.³

The inverse Gaussian distribution has several properties analogous to a Gaussian distribution. The name can be misleading: it is an inverse only in that, while the Gaussian describes a Brownian motion's level at a fixed time, the inverse Gaussian describes the distribution of the time a Brownian motion with positive drift takes to reach a fixed positive level. The relationship between the Gaussian and inverse Gaussian distributions is thus the same as the relationship between the binomial (number of successes for a fixed number of Bernoulli trials) and negative binomial (number of Bernoulli trials for a fixed number of successes) distributions.⁴

The y-axis reflections of the cumulant generating functions of the Gaussian and inverse Gaussian distributions are inverse of each other (i.e., the graphs of the two cumulant generating functions are reflections of each other across the line ⁠${\displaystyle y=-x}$⁠), a property that is also shared between the binomial and negative binomial distributions (after dividing their cumulant generating functions by their respective fixed parameter).⁴

To indicate that a random variable ⁠${\displaystyle X}$⁠ is inverse Gaussian-distributed with mean ⁠${\displaystyle \mu }$⁠ and shape parameter ⁠${\displaystyle \lambda }$⁠ we write ⁠${\displaystyle X\sim \operatorname {IG} (\mu ,\lambda )}$⁠.

The probability density function (pdf) of the inverse Gaussian distribution has a single parameter form given by

${\displaystyle f(x;\mu ,\mu ^{2})={\frac {\mu }{\sqrt {2\pi x^{3}}}}\exp {\biggl (}-{\frac {(x-\mu )^{2}}{2x}}{\biggr )}.}$

In this form, the mean and variance of the distribution are equal, ⁠${\displaystyle \mathbb {E} [X]=\operatorname {Var} (X)}$⁠.

Also, the cumulative distribution function (cdf) of the single parameter inverse Gaussian distribution is related to the standard normal distribution by

${\displaystyle {\begin{aligned}\Pr(X<x)&=\Phi (-z_{1})+e^{2\mu }\Phi (-z_{2}),\end{aligned}}}$

where ⁠${\displaystyle z_{1}={\frac {\mu }{x^{1/2}}}-x^{1/2}}$⁠, ⁠${\displaystyle z_{2}={\frac {\mu }{x^{1/2}}}+x^{1/2}}$⁠, and the ${\displaystyle \Phi }$ is the cdf of standard normal distribution. The variables ${\displaystyle z_{1}}$ and ${\displaystyle z_{2}}$ are related to each other by the identity ⁠${\displaystyle z_{2}^{2}=z_{1}^{2}+4\mu }$⁠.

In the single parameter form, the MGF simplifies to

${\displaystyle M(t)=\exp[\mu (1-{\sqrt {1-2t}})].}$

An inverse Gaussian distribution in double parameter form ${\displaystyle f(x;\mu ,\lambda )}$ can be transformed into a single parameter form ${\displaystyle f(y;\mu _{0},\mu _{0}^{2})}$ by appropriate scaling ⁠${\displaystyle y={\frac {\mu ^{2}x}{\lambda }}}$⁠, where ⁠${\displaystyle \mu _{0}=\mu ^{3}/\lambda }$⁠.

The above paragraph can be re-written as: if ⁠${\displaystyle Y=\lambda X/\mu ^{2}<}$⁠, then ⁠${\displaystyle Y\sim \operatorname {IG} (\lambda /\mu ,(\lambda /\mu )^{2})}$⁠⁵. This approach is better in the sense that it clearly shows dimensionless nature of the single parameter form (note that ⁠${\displaystyle \dim \lambda =\dim \mu =\dim x}$⁠). This property follows from a more general fact: if ${\displaystyle a>0}$ and ⁠${\displaystyle Y=aX}$⁠, then ⁠${\displaystyle Y\sim \operatorname {IG} (a\mu ,a\lambda )}$⁠².

The standard form of inverse Gaussian distribution is

${\displaystyle f(x;1,1)={\frac {1}{\sqrt {2\pi x^{3}}}}\exp {\biggl (}-{\frac {(x-1)^{2}}{2x}}{\biggr )}.}$

If ⁠${\displaystyle X_{i}}$⁠ has an ${\displaystyle \operatorname {IG} (\mu _{0}w_{i},\lambda _{0}w_{i}^{2})}$ distribution for ⁠${\displaystyle i=1,2,\dots ,n}$⁠ and all ⁠${\displaystyle X_{i}}$⁠ are independent, then

${\displaystyle S=\sum _{i=1}^{n}X_{i}\sim \operatorname {IG} \left(\mu _{0}\sum w_{i},\lambda _{0}\left(\sum w_{i}\right)^{2}\right).}$

The special case ${\displaystyle w_{i}=1/n}$ shows that the inverse Gaussian distribution is infinitely divisible.

Note that

${\displaystyle {\frac {\operatorname {Var} (X_{i})}{\operatorname {E} (X_{i})}}={\frac {\mu _{0}^{2}w_{i}^{2}}{\lambda _{0}w_{i}^{2}}}={\frac {\mu _{0}^{2}}{\lambda _{0}}}}$

is constant for all ⁠${\displaystyle i}$⁠. This is a necessary condition for the summation. Otherwise ⁠${\displaystyle S}$⁠ would not be Inverse Gaussian distributed.

For any ⁠${\displaystyle t>0}$⁠ it holds that

${\displaystyle X\sim \operatorname {IG} (\mu ,\lambda )\,\,\,\,\,\,\Rightarrow \,\,\,\,\,\,tX\sim \operatorname {IG} (t\mu ,t\lambda ).}$

The inverse Gaussian distribution is a two-parameter exponential family with natural parameters ⁠${\displaystyle -\lambda /(2\mu ^{2})}$⁠ and ⁠${\displaystyle -\lambda /2}$⁠, and natural statistics ⁠${\displaystyle X}$⁠ and ⁠${\displaystyle 1/X}$⁠.

For ${\displaystyle \lambda >0}$ fixed, it is also a single-parameter natural exponential family distribution⁶ where the base distribution has density

${\displaystyle h(x)={\sqrt {\frac {\lambda }{2\pi x^{3}}}}\exp \left(-{\frac {\lambda }{2x}}\right)\mathbb {1} _{[0,\infty )}(x)\,.}$

Indeed, with ⁠${\displaystyle \theta \leq 0}$⁠,

${\displaystyle p(x;\theta )={\frac {\exp(\theta x)h(x)}{\int \exp(\theta y)h(y)dy}}}$

is a density over the reals. Evaluating the integral, we get

${\displaystyle p(x;\theta )={\sqrt {\frac {\lambda }{2\pi x^{3}}}}\exp \left(-{\frac {\lambda }{2x}}+\theta x-{\sqrt {-2\lambda \theta }}\right)\mathbb {1} _{[0,\infty )}(x)\,.}$

Substituting ${\displaystyle \theta =-\lambda /(2\mu ^{2})}$ makes the above expression equal to ⁠${\displaystyle f(x;\mu ,\lambda )}$⁠.

Let the stochastic process ⁠${\displaystyle X_{t}}$⁠ be given by

${\displaystyle X_{0}=0\quad }$

${\displaystyle X_{t}=\nu t+\sigma W_{t}\quad \quad \quad \quad }$

where ⁠${\displaystyle W_{t}}$⁠ is a standard Brownian motion. That is, ⁠${\displaystyle X_{t}}$⁠ is a Brownian motion with drift ⁠${\displaystyle \nu >0}$⁠.

Then the first passage time for a fixed level ${\displaystyle \alpha >0}$ by ⁠${\displaystyle X_{t}}$⁠ is distributed according to an inverse-Gaussian:

${\displaystyle T_{\alpha }=\inf\{t>0\mid X_{t}=\alpha \}\sim \operatorname {IG} \left({\frac {\alpha }{\nu }},\left({\frac {\alpha }{\sigma }}\right)^{2}\right)={\frac {\alpha }{\sigma {\sqrt {2\pi x^{3}}}}}\exp {\biggl (}-{\frac {(\alpha -\nu x)^{2}}{2\sigma ^{2}x}}{\biggr )}}$

i.e

${\displaystyle P(T_{\alpha }\in (T,T+dT))={\frac {\alpha }{\sigma {\sqrt {2\pi T^{3}}}}}\exp {\biggl (}-{\frac {(\alpha -\nu T)^{2}}{2\sigma ^{2}T}}{\biggr )}dT}$

(cf. Schrödinger⁷ equation 19, Smoluchowski⁸, equation 8, and Folks⁵, equation 1).

A common special case of the above arises when the Brownian motion has no drift. In that case, parameter ⁠${\displaystyle \mu }$⁠ tends to infinity, and the first passage time for fixed level ⁠${\displaystyle \alpha }$⁠ has probability density function

${\displaystyle f\left(x;0,\left({\frac {\alpha }{\sigma }}\right)^{2}\right)={\frac {\alpha }{\sigma {\sqrt {2\pi x^{3}}}}}\exp \left(-{\frac {\alpha ^{2}}{2\sigma ^{2}x}}\right)}$

(see also Bachelier^{9: 74 10: 39}). This is a Lévy distribution with parameters ${\displaystyle c=\left({\frac {\alpha }{\sigma }}\right)^{2}}$ and ⁠${\displaystyle \mu =0}$⁠.

The model where

${\displaystyle X_{i}\sim \operatorname {IG} (\mu ,\lambda w_{i}),\,\,\,\,\,\,i=1,2,\ldots ,n}$

with all ⁠${\displaystyle w_{i}}$⁠ known, ⁠${\displaystyle (\mu ,\lambda )}$⁠ unknown and all ⁠${\displaystyle X_{i}}$⁠ independent has the following likelihood function:

${\displaystyle L(\mu ,\lambda )=\left({\frac {\lambda }{2\pi }}\right)^{\frac {n}{2}}\left(\prod _{i=1}^{n}{\frac {w_{i}}{X_{i}^{3}}}\right)^{\frac {1}{2}}\exp \left({\frac {\lambda }{\mu }}\sum _{i=1}^{n}w_{i}-{\frac {\lambda }{2\mu ^{2}}}\sum _{i=1}^{n}w_{i}X_{i}-{\frac {\lambda }{2}}\sum _{i=1}^{n}w_{i}{\frac {1}{X_{i}}}\right).}$

Solving the likelihood equation yields the following maximum likelihood estimates

${\displaystyle {\widehat {\mu }}={\frac {\sum _{i=1}^{n}w_{i}X_{i}}{\sum _{i=1}^{n}w_{i}}},\,\,\,\,\,\,\,\,{\frac {1}{\widehat {\lambda }}}={\frac {1}{n}}\sum _{i=1}^{n}w_{i}\left({\frac {1}{X_{i}}}-{\frac {1}{\widehat {\mu }}}\right).}$

${\displaystyle {\widehat {\mu }}}$ and ${\displaystyle {\widehat {\lambda }}}$ are independent and

${\displaystyle {\widehat {\mu }}\sim \operatorname {IG} \left(\mu ,\lambda \sum _{i=1}^{n}w_{i}\right),\qquad {\frac {n}{\widehat {\lambda }}}\sim {\frac {1}{\lambda }}\chi _{n-1}^{2}.}$

The following algorithm may be used.¹¹

Generate a random variate from a normal distribution with mean ⁠${\displaystyle 0}$⁠ and standard deviation equal ⁠${\displaystyle 1}$⁠

${\displaystyle \displaystyle \nu \sim N(0,1).}$

Square the value

${\displaystyle \displaystyle y=\nu ^{2}}$

and use the relation

${\displaystyle x=\mu +{\frac {\mu ^{2}y}{2\lambda }}-{\frac {\mu }{2\lambda }}{\sqrt {4\mu \lambda y+\mu ^{2}y^{2}}}.}$

Generate another random variate, this time sampled from a uniform distribution between ⁠${\displaystyle 0}$⁠ and ⁠${\displaystyle 1}$⁠

${\displaystyle \displaystyle z\sim U(0,1).}$

If ${\displaystyle z\leq {\frac {\mu }{\mu +x}}}$ then return ${\displaystyle \displaystyle x}$ else return ${\displaystyle {\frac {\mu ^{2}}{x}}.}$

Sample code in Java:

public double inverseGaussian(double mu, double lambda) {
    Random rand = new Random();
    double v = rand.nextGaussian();  // Sample from a normal distribution with a mean of 0 and 1 standard deviation
    double y = v * v;
    double x = mu + (mu * mu * y) / (2 * lambda) - (mu / (2 * lambda)) * Math.sqrt(4 * mu * lambda * y + mu * mu * y * y);
    double test = rand.nextDouble();  // Sample from a uniform distribution between 0 and 1
    if (test <= (mu) / (mu + x))
        return x;
    else
        return (mu * mu) / x;
}

And to plot Wald distribution in Python using matplotlib and NumPy:

import matplotlib.pyplot as plt
import numpy as np

h = plt.hist(np.random.wald(3, 2, 100000), bins = 200, density = True)

plt.show()

• If ⁠${\displaystyle X\sim \operatorname {IG} (\mu ,\lambda )}$⁠, then ⁠${\displaystyle kX\sim \operatorname {IG} (k\mu ,k\lambda )}$⁠ for any number ⁠${\displaystyle k>0}$⁠.¹
• If ${\displaystyle X_{i}\sim \operatorname {IG} (\mu ,\lambda )}$ then ⁠${\displaystyle {1}}$⁠.
• If ${\displaystyle X_{i}\sim \operatorname {IG} (\mu ,\lambda )}$ for ${\displaystyle i=1,\ldots ,n}$ then ⁠${\displaystyle {\bar {X}}\sim \operatorname {IG} (\mu ,n\lambda )}$⁠.
• If ${\displaystyle X_{i}\sim \operatorname {IG} (\mu _{i},2\mu _{i}^{2})}$ then ⁠${\displaystyle \textstyle \sum _{i=1}^{n}X_{i}\sim \operatorname {IG} \left(\sum _{i=1}^{n}\mu _{i},2\left(\sum _{i=1}^{n}\mu _{i}\right)^{2}\right)}$⁠.
• If ⁠${\displaystyle X\sim \operatorname {IG} (\mu ,\lambda )}$⁠, then ⁠${\displaystyle \lambda (X-\mu )^{2}/\mu ^{2}X\sim \chi ^{2}(1)}$⁠.¹²

The convolution of an inverse Gaussian distribution (a Wald distribution) and an exponential (an ex-Wald distribution) is used as a model for response times in psychology,¹³ with visual search as one example.¹⁴

This distribution appears to have been first derived in 1900 by Louis Bachelier⁹¹⁰ as the time a stock reaches a certain price for the first time. In 1915 it was used independently by Erwin Schrödinger⁷ and Marian v. Smoluchowski⁸ as the time to first passage of a Brownian motion. In the field of reproduction modeling it is known as the Hadwiger function, after Hugo Hadwiger who described it in 1940.¹⁵ Abraham Wald re-derived this distribution in 1944¹⁶ as the limiting form of a sample in a sequential probability ratio test. The name inverse Gaussian was proposed by Maurice Tweedie in 1945.⁴ Tweedie investigated this distribution in 1956¹⁷ and 1957²³ and established some of its statistical properties. The distribution was extensively reviewed by Folks and Chhikara in 1978.⁵

Assuming that the time intervals between occurrences of a random phenomenon follow an inverse Gaussian distribution, the probability distribution for the number of occurrences of this event within a specified time window is referred to as rated inverse Gaussian.¹⁸ While, first and second moment of this distribution are calculated, the derivation of the moment generating function remains an open problem.

Despite the simple formula for the probability density function, numerical probability calculations for the inverse Gaussian distribution nevertheless require special care to achieve full machine accuracy in floating point arithmetic for all parameter values.¹⁹ Functions for the inverse Gaussian distribution are provided for the R programming language by several packages including rmutil,²⁰²¹ SuppDists,²² STAR,²³ invGauss,²⁴ LaplacesDemon,²⁵ and statmod.²⁶