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 (0,∞).
Its probability density function is given by
f(x;\mu,λ)=\sqrt
λ | \expl(- | |
2\pix3 |
λ(x-\mu)2 | |
2\mu2x |
r)
for x > 0, where
\mu>0
λ>0
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.
Its cumulant generating function (logarithm of the characteristic function) is the inverse of the cumulant generating function of a Gaussian random variable.
To indicate that a random variable X is inverse Gaussian-distributed with mean μ and shape parameter λ we write
X\sim\operatorname{IG}(\mu,λ)
The probability density function (pdf) of the inverse Gaussian distribution has a single parameter form given by
f(x;\mu,\mu2) =
\mu | |
\sqrt{2\pix3 |
In this form, the mean and variance of the distribution are equal,
E[X]=Var(X).
Also, the cumulative distribution function (cdf) of the single parameter inverse Gaussian distribution is related to the standard normal distribution by
\begin{align} \Pr(X<x)&=\Phi(-z1)+e2\Phi(-z2), \end{align}
where
z1=
\mu | |
x1/2 |
-x1/2
z2=
\mu | |
x1/2 |
+x1/2,
\Phi
z1
z2
2 | |
z | |
2 |
=
2 | |
z | |
1 |
+4\mu.
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,λ)
f(y;\mu0,\mu
2) | |
0 |
y=
\mu2x | |
λ |
,
\mu0=\mu3/λ.
The standard form of inverse Gaussian distribution is
f(x;1,1) =
1 | |
\sqrt{2\pix3 |
If Xi has an
\operatorname{IG}(\mu0wi,λ0
2 | |
w | |
i |
)
n | |
S=\sum | |
i=1 |
Xi \sim \operatorname{IG}\left(\mu0\sumwi,λ0\left(\sumwi\right)2\right).
Note that
\operatorname{Var | |
(X |
i)}{\operatorname{E}(Xi)}=
| = | ||||||||||||||||
|
| |||||||
λ0 |
is constant for all i. This is a necessary condition for the summation. Otherwise S would not be Inverse Gaussian distributed.
For any t > 0 it holds that
X\sim\operatorname{IG}(\mu,λ) ⇒ tX\sim\operatorname{IG}(t\mu,tλ).
The inverse Gaussian distribution is a two-parameter exponential family with natural parameters −λ/(2μ2) and −λ/2, and natural statistics X and 1/X.
For
λ>0
h(x)=\sqrt{
λ | |
2\pix3 |
}\exp\left(-
λ | |
2x |
\right)1[0,infty)(x).
Indeed, with
\theta\le0
p(x;\theta)=
\exp(\thetax)h(x) | |
\int\exp(\thetay)h(y)dy |
is a density over the reals. Evaluating the integral, we get
p(x;\theta)=\sqrt{
λ | |
2\pix3 |
}\exp\left(-
λ | |
2x |
+\thetax-\sqrt{-2λ\theta}\right)1[0,infty)(x).
Substituting
\theta=-λ/(2\mu2)
f(x;\mu,λ)
Let the stochastic process Xt be given by
X0=0
Xt=\nut+\sigmaWt
where Wt is a standard Brownian motion. That is, Xt is a Brownian motion with drift
\nu>0
Then the first passage time for a fixed level
\alpha>0
T\alpha=inf\{t>0\midXt=\alpha\}\sim\operatorname{IG}\left(
\alpha\nu, | |||
|
\right)2\right) =
\alpha | |
\sigma\sqrt{2\pix3 |
P(T\alpha\in(T,T+dT))=
\alpha | |
\sigma\sqrt{2\piT3 |
(cf. Schrödinger equation 19, Smoluchowski, equation 8, and Folks, equation 1).
Suppose that we have a Brownian motion
Xt
\nu
Xt=\nut+\sigmaWt, X(0)=x0
And suppose that we wish to find the probability density function for the time when the process first hits some barrier
\alpha>x0
p(t,x)
{\partialp\over{\partialt}}+\nu{\partialp\over{\partialx}}={1\over{2}}\sigma2{\partial2p\over{\partialx2
where
\delta( ⋅ )
p(t,\alpha)=0
\varphi(t,x)
\varphi(t,x)={1\over{\sqrt{2\pi\sigma2t}}}\exp\left[-{(x-x0-\nut)2\over{2\sigma2t}}\right]
Define a point
m
m>\alpha
p(0,x)=\delta(x-x0)-A\delta(x-m)
where
A
p(t,x)={1\over{\sqrt{2\pi\sigma2t}}}\left\{\exp\left[-{(x-x0-\nut)2\over{2\sigma2t}}\right]-A\exp\left[-{(x-m-\nut)2\over{2\sigma2t}}\right]\right\}
Now we must determine the value of
A
(\alpha-x0-\nut)2=-2\sigma2tlogA+(\alpha-m-\nut)2
At
p(0,\alpha)
(\alpha-x0)2=(\alpha-m)2\impliesm=2\alpha-x0
A=
2\nu(\alpha-x0)/\sigma2 | |
e |
Therefore, the full solution to the BVP is:
p(t,x)={1\over{\sqrt{2\pi\sigma2t}}}\left\{\exp\left[-{(x-x0-\nut)2\over{2\sigma2t}}\right]-
2\nu(\alpha-x0)/\sigma2 | |
e |
\exp\left[-{(x+x0-2\alpha-\nut)2\over{2\sigma2t}}\right]\right\}
Now that we have the full probability density function, we are ready to find the first passage time distribution
f(t)
S(t)
\begin{aligned}S(t)&=
\alpha | |
\int | |
-infty |
p(t,x)dx\ &=\Phi\left({\alpha-x0-\nut\over{\sigma\sqrt{t}}}\right)-
2\nu(\alpha-x0)/\sigma2 | |
e |
\Phi\left({-\alpha+x0-\nut\over{\sigma\sqrt{t}}}\right)\end{aligned}
where
\Phi( ⋅ )
\alpha
t
f(t)
\begin{aligned}f(t)&=-{dS\over{dt}}\ &={(\alpha-x0)\over{\sqrt{2\pi\sigma2t3
Assuming that
x0=0
f(t)={\alpha\over{\sqrt{2\pi\sigma2t3
A common special case of the above arises when the Brownian motion has no drift. In that case, parameter μ tends to infinity, and the first passage time for fixed level α has probability density function
f\left(x;0,\left(
\alpha | |
\sigma |
\right)2\right) =
\alpha | |
\sigma\sqrt{2\pix3 |
(see also Bachelier). This is a Lévy distribution with parameters
c=\left( | \alpha |
\sigma |
\right)2
\mu=0
The model where
Xi\sim\operatorname{IG}(\mu,λwi),i=1,2,\ldots,n
with all wi known, (μ, λ) unknown and all Xi independent has the following likelihood function
L(\mu,λ)= \left(
λ | |
2\pi |
| ||||
\right) |
n | |
\prod | |
i=1 |
wi | ||||||
|
| |||||
\right) | \exp\left( |
λ | |
\mu |
n | |
\sum | |
i=1 |
wi-
λ | |
2\mu2 |
n | |
\sum | |
i=1 |
wiXi-
λ | |
2 |
n | |
\sum | |
i=1 |
wi
1{X | |
i} |
\right).
Solving the likelihood equation yields the following maximum likelihood estimates
\widehat{\mu}=
| ||||||||||
|
,
1 | |
\widehat{λ |
\widehat{\mu}
\widehat{λ}
\widehat{\mu}\sim\operatorname{IG}\left(\mu,λ
n | |
\sum | |
i=1 |
wi\right),
n | |
\widehat{λ |
The following algorithm may be used.
Generate a random variate from a normal distribution with mean 0 and standard deviation equal 1Sample code in Java:\displaystyle\nu\simN(0,1).
Square the value
\displaystyley=\nu2
and use the relation
x=\mu+
\mu2y 2λ -
\mu 2λ \sqrt{4\muλy+\mu2y2}.
Generate another random variate, this time sampled from a uniform distribution between 0 and 1
\displaystylez\simU(0,1).
If
then returnz\le
\mu \mu+x else return\displaystyle x
\mu2 x .
And to plot Wald distribution in Python using matplotlib and NumPy:
h = plt.hist(np.random.wald(3, 2, 100000), bins=200, density=True)
plt.show
X\sim\operatorname{IG}(\mu,λ)
kX\sim\operatorname{IG}(k\mu,kλ)
k>0.
Xi\sim\operatorname{IG}(\mu,λ)
n | |
\sum | |
i=1 |
Xi\sim\operatorname{IG}(n\mu,n2λ)
Xi\sim\operatorname{IG}(\mu,λ)
i=1,\ldots,n
\bar{X}\sim\operatorname{IG}(\mu,nλ)
Xi\sim\operatorname{IG}(\mui,2
2 | |
\mu | |
i) |
n | |
\sum | |
i=1 |
Xi\sim
n | |
\operatorname{IG}\left(\sum | |
i=1 |
\mui,2\left(
n | |
\sum | |
i=1 |
\mui\right)2\right)
X\sim\operatorname{IG}(\mu,λ)
λ(X-\mu)2/\mu2X\sim\chi2(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.[2]
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.[3] 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.[4] Tweedie investigated this distribution in 1956[5] and 1957[6] [7] and established some of its statistical properties. The distribution was extensively reviewed by Folks and Chhikara in 1978.
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.[8] Functions for the inverse Gaussian distribution are provided for the R programming language by several packages including rmutil,[9] [10] SuppDists,[11] STAR,[12] invGauss,[13] LaplacesDemon,[14] and statmod.[15]