Inverse Gaussian distribution

Let the stochastic process X_t be given by

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

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

where W_t is a standard Brownian motion. That is, X_t is a Brownian motion with drift ${\displaystyle \nu >0}$.

Then the first passage time for a fixed level ${\displaystyle \alpha >0}$ by 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 )}}$

(cf. Schrödinger[2] equation 19, Smoluchowski[3], equation 8, and Folks[4], equation 1).

The model where

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

with all w_i known, (μ, λ) unknown and all 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.[7]

Generate a random variate from a normal distribution with mean 0 and standard deviation equal 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 0 and 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;
}

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

• If ${\displaystyle X\sim \operatorname {IG} (\mu ,\lambda )}$, then ${\displaystyle kX\sim \operatorname {IG} (k\mu ,k\lambda )}$ for any number ${\displaystyle k>0.}$[1]
• If ${\displaystyle X_{i}\sim \operatorname {IG} (\mu ,\lambda )\,}$ then ${\displaystyle \sum _{i=1}^{n}X_{i}\sim \operatorname {IG} (n\mu ,n^{2}\lambda )\,}$
• 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 \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)\,}$

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,[8] with visual search as one example.[9]

This distribution appears to have been first derived in 1900 by Louis Bachelier[5][6] as the time a stock reaches a certain price for the first time. In 1915 it was used independently by Erwin Schrödinger[2] and Marian v. Smoluchowski[3] 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.[10] Abraham Wald re-derived this distribution in 1944[11] as the limiting form of a sample in a sequential probability ratio test. The name inverse Gaussian was proposed by Maurice Tweedie in 1945.[12] Tweedie investigated this distribution in 1956[13] and 1957[14] [15] and established some of its statistical properties. The distribution was extensively reviewed by Folks and Chhikara in 1978.[4]

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.[16] Functions for the inverse Gaussian distribution are provided for the R programming language by several packages including rmutil,[17][18] SuppDists,[19] STAR,[20] invGauss,[21] LaplacesDemon,[22] and statmod.[23]

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

• Høyland, Arnljot; Rausand, Marvin (1994). System Reliability Theory. New York: Wiley. ISBN 978-0-471-59397-3.
• Seshadri, V. (1993). The Inverse Gaussian Distribution. Oxford University Press. ISBN 978-0-19-852243-0.

• Inverse Gaussian Distribution in Wolfram website.