Inverse Gaussian distribution

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 Template:Tmath.

Its probability density function is given by

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

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

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.Template:R

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 Template:Tmath), a property that is also shared between the binomial and negative binomial distributions (after dividing their cumulant generating functions by their respective fixed parameter).Template:R

To indicate that a random variable Template:Tmath is inverse Gaussian-distributed with mean Template:Tmath and shape parameter Template:Tmath we write Template:Tmath.

Properties

Single parameter form

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

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

In this form, the mean and variance of the distribution are equal, Template:Tmath.

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

$\begin{array}{cc}Pr(X<x)& =\mathrm{\Phi }(-z_1)+e^{2\mu }\mathrm{\Phi }(-z_2),\end{array}$

where Template:Tmath, Template:Tmath, and the $\mathrm{\Phi }$ is the cdf of standard normal distribution. The variables $z_1$ and $z_2$ are related to each other by the identity Template:Tmath.

In the single parameter form, the MGF simplifies to

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

An inverse Gaussian distribution in double parameter form $f(x;\mu ,\lambda )$ can be transformed into a single parameter form $f(y;{\mu }_0,{\mu }_0^2)$ by appropriate scaling Template:Tmath, where Template:Tmath.

The above paragraph can be re-written as: if Template:Tmath, then Template:TmathTemplate:R. This approach is better in the sense that it clearly shows dimensionless nature of the single parameter form (note that Template:Tmath). This property follows from a more general fact: if $a>0$ and Template:Tmath, then Template:TmathTemplate:R.

The standard form of inverse Gaussian distribution is

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

Summation

If Template:Tmath has an $\mathrm{IG}({\mu }_0{w}_i,{\lambda }_0{w}_i^2)$ distribution for Template:Tmath and all Template:Tmath are independent, then

$S=\sum _{i=1}^n{X}_i\sim \mathrm{IG}({\mu }_0\sum w_i,{\lambda }_0{(\sum w_i)}^2).$

The special case $w_i=1/n$ shows that the inverse Gaussian distribution is infinitely divisible.

Note that

$\frac{\mathrm{Var}(X_i)}{\mathrm{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 Template:Tmath. This is a necessary condition for the summation. Otherwise Template:Tmath would not be Inverse Gaussian distributed.

Scaling

For any Template:Tmath it holds that

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

Exponential family

The inverse Gaussian distribution is a two-parameter exponential family with natural parameters Template:Tmath and Template:Tmath, and natural statistics Template:Tmath and Template:Tmath.

For $\lambda >0$ fixed, it is also a single-parameter natural exponential family distributionTemplate:R where the base distribution has density

$h(x)=\sqrt{\frac{\lambda }{2\pi x^3}}\exp (-\frac{\lambda }{2x}){\mathrm{𝟙}}_{[0,\mathrm{\infty })}(x).$

Indeed, with Template:Tmath,

$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

$p(x;\theta )=\sqrt{\frac{\lambda }{2\pi x^3}}\exp (-\frac{\lambda }{2x}+\theta x-\sqrt{-2\lambda \theta }){\mathrm{𝟙}}_{[0,\mathrm{\infty })}(x).$

Substituting $\theta =-\lambda /(2{\mu }^2)$ makes the above expression equal to Template:Tmath.

Relationship with Brownian motion

Example of stopped random walks with Template:Tmath. The upper figure shows the histogram of waiting times, along with the prediction according to inverse gaussian distribution. The lower figure shows the trajectories.

Let the stochastic process Template:Tmath be given by

$X_0=0$

$X_t=\nu t+\sigma W_t$

where Template:Tmath is a standard Brownian motion. That is, Template:Tmath is a Brownian motion with drift Template:Tmath.

Then the first passage time for a fixed level $\alpha >0$ by Template:Tmath is distributed according to an inverse-Gaussian:

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

i.e

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

(cf. SchrödingerTemplate:R equation 19, SmoluchowskiTemplate:R, equation 8, and FolksTemplate:R, equation 1).

Template:Collapse top Suppose that we have a Brownian motion $X_t$ with drift $\nu$ defined by:

$X_t=\nu t+\sigma W_t,X(0)=x_0$

And suppose that we wish to find the probability density function for the time when the process first hits some barrier $\alpha >x_0$ - known as the first passage time. The Fokker–Planck equation describing the evolution of the probability distribution $p(t,x)$ is:

$\frac{\partial p}{\partial t}+\nu \frac{\partial p}{\partial x}=\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}p}{\partial x^2},\{\begin{array}{cc}p(0,x)& =\delta (x-x_0)\\ p(t,\alpha )& =0\end{array}$

where $\delta (\cdot )$ is the Dirac delta function. This is a boundary value problem (BVP) with a single absorbing boundary condition $p(t,\alpha )=0$, which may be solved using the method of images. Based on the initial condition, the fundamental solution to the Fokker–Planck equation, denoted by $\phi (t,x)$, is:

$\phi (t,x)=\frac{1}{\sqrt{2\pi {\sigma }^{2}t}}\exp [-\frac{(x-x_0-\nu t{)}^2}{2{\sigma }^{2}t}]$

Define a point Template:Tmath, such that Template:Tmath. This will allow the original and mirror solutions to cancel out exactly at the barrier at each instant in time. This implies that the initial condition should be augmented to become:

$p(0,x)=\delta (x-x_0)-A\delta (x-m)$

where $A$ is a constant. Due to the linearity of the BVP, the solution to the Fokker–Planck equation with this initial condition is:

$p(t,x)=\frac{1}{\sqrt{2\pi {\sigma }^{2}t}}\{\exp [-\frac{(x-x_0-\nu t{)}^2}{2{\sigma }^{2}t}]-A\exp [-\frac{(x-m-\nu t{)}^2}{2{\sigma }^{2}t}]\}$

Now we must determine the value of Template:Tmath. The fully absorbing boundary condition implies that:

$(\alpha -x_0-\nu t{)}^2=-2{\sigma }^{2}t\log A+(\alpha -m-\nu t{)}^2$

At Template:Tmath, we have that Template:Tmath. Substituting this back into the above equation, we find that:

$A=e^{2\nu (\alpha -x_0)/{\sigma }^2}$

Therefore, the full solution to the BVP is:

$p(t,x)=\frac{1}{\sqrt{2\pi {\sigma }^{2}t}}\{\exp [-\frac{(x-x_0-\nu t{)}^2}{2{\sigma }^{2}t}]-e^{2\nu (\alpha -x_0)/{\sigma }^2}\exp [-\frac{(x+x_0-2\alpha -\nu t{)}^2}{2{\sigma }^{2}t}]\}$

Now that we have the full probability density function, we are ready to find the first passage time distribution Template:Tmath. The simplest route is to first compute the survival function Template:Tmath, which is defined as:

$\begin{array}{cc}S(t)& ={\int }_{-\mathrm{\infty }}^{\alpha }p(t,x)dx\\ & =\mathrm{\Phi }\left(\frac{\alpha -x_0-\nu t}{\sigma \sqrt{t}}\right)-e^{2\nu (\alpha -x_0)/{\sigma }^2}\mathrm{\Phi }\left(\frac{-\alpha +x_0-\nu t}{\sigma \sqrt{t}}\right)\end{array}$

where $\mathrm{\Phi }(\cdot )$ is the cumulative distribution function of the standard normal distribution. The survival function gives us the probability that the Brownian motion process has not crossed the barrier $\alpha$ at some time Template:Tmath. Finally, the first passage time distribution $f(t)$ is obtained from the identity:

$\begin{array}{cc}f(t)& =-\frac{dS}{dt}\\ & =\frac{(\alpha -x_0)}{\sqrt{2\pi {\sigma }^2{t}^3}}{e}^{-(\alpha -x_0-\nu t{)}^2/2{\sigma }^{2}t}\end{array}$

Assuming that $x_0=0$, the first passage time follows an inverse Gaussian distribution:

$f(t)=\frac{\alpha }{\sqrt{2\pi {\sigma }^2{t}^3}}{e}^{-(\alpha -\nu t{)}^2/2{\sigma }^{2}t}\sim \text{IG}[\frac{\alpha }{\nu },{\left(\frac{\alpha }{\sigma }\right)}^2]$

When drift is zero

A common special case of the above arises when the Brownian motion has no drift. In that case, parameter Template:Tmath tends to infinity, and the first passage time for fixed level Template:Tmath has probability density function

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

(see also BachelierTemplate:RTemplate:R). This is a Lévy distribution with parameters $c={\left(\frac{\alpha }{\sigma }\right)}^2$ and Template:Tmath.

Maximum likelihood

The model where

$X_i\sim \mathrm{IG}(\mu ,\lambda w_i),i=1,2,\dots ,n$

with all Template:Tmath known, Template:Tmath unknown and all Template:Tmath independent has the following likelihood function:

$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 (\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}).$

Solving the likelihood equation yields the following maximum likelihood estimates

$\hat{\mu }=\frac{\sum _{i=1}^n{w}_i{X}_i}{\sum _{i=1}^n{w}_i},\frac{1}{\hat{\lambda }}=\frac{1}{n}\sum _{i=1}^n{w}_i(\frac{1}{X_i}-\frac{1}{\hat{\mu }}).$

$\hat{\mu }$ and $\hat{\lambda }$ are independent and

$\hat{\mu }\sim \mathrm{IG}(\mu ,\lambda \sum _{i=1}^n{w}_i),\frac{n}{\hat{\lambda }}\sim \frac{1}{\lambda }{\chi }_{n-1}^2.$

Sampling from an inverse-Gaussian distribution

The following algorithm may be used.Template:R

Generate a random variate from a normal distribution with mean Template:Tmath and standard deviation equal Template:Tmath

$\nu \sim N(0,1).$

Square the value

$y={\nu }^2$

and use the relation

$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 Template:Tmath and Template:Tmath

$z\sim U(0,1).$

If $z\le \frac{\mu }{\mu +x}$ then return $x$ else return $\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;
}

Wald distribution using Python with aid of matplotlib and NumPy

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()

Related distributions

• If Template:Tmath, then Template:Tmath for any number Template:Tmath.Template:R
• If $X_i\sim \mathrm{IG}(\mu ,\lambda )$ then Template:Tmath.
• If $X_i\sim \mathrm{IG}(\mu ,\lambda )$ for $i=1,\dots ,n$ then Template:Tmath.
• If $X_i\sim \mathrm{IG}({\mu }_i,2{\mu }_i^2)$ then Template:Tmath.
• If Template:Tmath, then Template:Tmath.

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,Template:R with visual search as one example.

History

This distribution appears to have been first derived in 1900 by Louis BachelierTemplate:R as the time a stock reaches a certain price for the first time. In 1915 it was used independently by Erwin SchrödingerTemplate:R and Marian v. SmoluchowskiTemplate:R 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 1944Template:R 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.Template:R

Rated inverse Gaussian distribution

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.

Numeric computation and software

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.

See also

• Generalized inverse Gaussian distribution
• Tweedie distribution – The inverse Gaussian distribution is a member of the family of Tweedie exponential dispersion models
• Stopping time

Further reading

• Template:Cite book
• Template:Cite book

External links

• Inverse Gaussian Distribution in Wolfram website.