# Can't compute CDF for Inverse Gaussian distribution

I am trying to implement in Python the CDF of the Inverse Gaussian distribution:

Inverse Gaussian pdf: $$f(x) = \sqrt{\frac{\lambda}{2\pi x^3}}e^{-\frac{\lambda(x-\mu)^2}{2\mu^2x}}$$

Inverse Gaussian cdf: $$F(x) = \Phi\left(\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu}-1\right)\right) + e^{\frac{2\lambda}{\mu}}\Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu}+1\right)\right)$$

Where $$\Phi$$ is the standard Normal distribution cdf

import numpy as np
from scipy.stats import norm

def igxy(mu,lamb):
def _ig(x,mu,lamb):
arg = lamb / (2 * np.pi * x ** 3)
return np.sqrt(arg) * np.exp((-lamb * (x - mu) ** 2) / (2 * mu ** 2 * x))
x = np.linspace(0.0001,5,2000)
y = _ig(x,mu,lamb)
return x,y

def igcdf(x: float, mu: float, lamb:float):
return np.nan_to_num(
norm.cdf(np.sqrt(lamb/x) * (x/mu -1)) * np.exp(2*lamb/mu) * norm.cdf(-np.sqrt(lamb/x)*(x/mu + 1))
)

mu = 1.1
lamb = 3
x,y = igxy(mu,lamb)
plt.figure(figsize=(10,5),dpi = 150)
delta = (x[-1] - x[0]) / len(x)
plt.plot(x, y, 'k-', lw=2)
plt.plot(x,np.cumsum(y)*delta)
plt.plot(x,[igcdf(el,mu,lamb) for el in x])

plt.legend(["Inverse Gaussian PDF","Expected Inverse Gaussian CDF","Inverse Gaussian CDF"])


Now, the first problem I have is that the distribution I obtain in code doesn't seem to be a cdf to me (by looking at the graph). Moreover I need to work with values of $$\lambda\sim1500$$, and in Python this translates in overflow errors when it has to compute the exponential, how can I manage this kind of situations?

Thanks

1. Finally in order to avoid overflows I found the following solution:
import numpy as np
import matplotlib.pyplot as plt
from typing import Union
from scipy.stats import invgauss,norm

_MAX_EXPONENT = 700
def _exp(n):
return np.exp(min(n,_MAX_EXPONENT))

def igpdf(x: np.array , mu: float, lamb: float) -> np.array:
arg = lamb / (2 * np.pi * x ** 3)
return np.sqrt(arg) * np.exp((-lamb * (x - mu) ** 2) / (2 * mu ** 2 * x))

def igcdf(x: np.array, mu: float, lamb:float):
return norm.cdf(np.sqrt(lamb/x) * (x/mu -1)) + _exp(2 * lamb / mu) * norm.cdf(-np.sqrt(lamb / x) * (x / mu + 1))

mu = 1.1
lamb = 1700
x = np.linspace(0.001,2.000,1000)

plt.figure(figsize=(10,5),dpi = 150)
plt.plot(x, igpdf(x,mu,lamb), 'k-', lw=2)
plt.plot(x, igcdf(x,mu,lamb), 'b-', lw=2)

plt.legend(["Inverse Gaussian PDF","Inverse Gaussian CDF"])


Which remains a better solution w.r.t. to using scipy:

from scipy.stats import invgauss
plt.figure(figsize=(10,5),dpi = 150)
plt.plot(x,invgauss.pdf(x,mu = mu/lamb,scale = lamb))
plt.plot(x,invgauss.cdf(x,mu = mu/lamb,scale = lamb))
plt.legend(["Scipy IG PDF","Scipy IG CDF"])

/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/scipy/stats/_continuous_distns.py:3613: RuntimeWarning: overflow encountered in exp
C1 += np.exp(1.0/mu) * _norm_cdf(-fac*(x+mu)/mu) * np.exp(1.0/mu)
/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/scipy/stats/_continuous_distns.py:3613: RuntimeWarning: invalid value encountered in multiply
C1 += np.exp(1.0/mu) * _norm_cdf(-fac*(x+mu)/mu) * np.exp(1.0/mu)


Now, the problem that I have is that I'm trying to parametrise the Inverse Gaussian distribution in order to train a model through Maximum Likelihood Estimation, to do so I had to compute the gradient and hessian matrix of the following function w.r.t each parameter.

$$\mathcal{L}(\mu,\lambda)=-\eta_n\log(1-P(t))$$

For example the second order derivative w.r.t $$\mu$$ looks like this:

$$\frac{\partial^2}{\partial^2\mu}\mathcal{L}(\mu,\lambda) = \\ \frac{-k \eta \left(e^{-\frac{k (t-\mu )^2}{2 \mu ^2 t}} \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu -t)}{\mu }\right)+ e^{\frac{2 k}{\mu }} \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu +t)}{\mu }\right)-1\right)\right) \left(4 (k+\mu ) e^{\frac{k (\mu +t)^2}{2 \mu ^2 t}} \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu +t)}{\mu }\right)-1\right)+ \sqrt{\frac{2}{\pi }} t \sqrt{\frac{k}{t}}\right)-4 k e^{\frac{4 k}{\mu }} \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu +t)}{\mu }\right)-1\right)^2\right)} {\mu ^4 \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu -t)}{\mu }\right)+ e^{\frac{2 k}{\mu }} \left(\Phi \left(\frac{\sqrt{\frac{k}{t}} (\mu +t)}{\mu }\right)-1\right)\right)^2}$$

def _right_censoring_derivative_thetatheta(k, t, mu, etan, xt):
numerator_f1 = k*etan
numerator_f2_add1_f1 = -_exp(-(k * (t - mu) ** 2) / (2 * t * mu ** 2))
numerator_f2_add1_f2 = _norm_cdf((np.sqrt(k/t)*(-t+mu))/mu) + _exp(2 * k / mu) * (-1 + _norm_cdf(np.sqrt(k / t) * (t + mu) / mu))
numerator_f2_add1_f3 = np.sqrt(2/np.pi) * np.sqrt(k/t) * t + 4 * _exp(k * (t + mu) ** 2 / (2 * t * mu ** 2)) * (k + mu) * (-1 + _norm_cdf((np.sqrt(k / t) * (t + mu)) / mu))
numerator_f2_add2 = 4 * _exp(4 * k / mu) * k * (-1 + _norm_cdf((np.sqrt(k / t) * (t + mu)) / mu)) ** 2
denominator = mu ** 4 * (_norm_cdf(np.sqrt(k/t)*(-t+mu)/mu) + _exp(2 * k / mu) * (-1 + _norm_cdf(np.sqrt(k / t) * (t + mu) / mu))) ** 2
xtx = np.dot(xt.T,xt)
return numerator/denominator * xtx


And I'm finding a lot of overflow errors when actually computing these values, it would be a good idea to "clip" the overflow values between an upper bound(e.g. 1e100) and a lower bound (1e-100)? Or are there smarter ways to overcome this kind of issues?

Thanks again,

• (x/mu -1)) * np.exp(2*lamb/mu should be a +, no? Commented Nov 22, 2020 at 13:36
• what's the order of magnitude of $\mu$? Commented Nov 22, 2020 at 13:39
• @Cam.Davidson.Pilon You're totally right, I have no idea how I missed that. $\mu$ is always $0.6<\mu<1.2$ Commented Nov 22, 2020 at 13:51

1. We fixed the first question by finding a transcription error.

2. The second question about overflows can be solved as follows. The terms:

$$e^{\frac{2\lambda}{\mu}}\Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu}+1\right)\right)$$

are the multiplication of a very large number and a very small number. In fact, the second term, for the values you provided, evaluates to 0 on my machine:

lambda_ = 1500
mu = 1
x = 1
stats.norm.cdf(-np.sqrt(lambda_/mu) * (x/mu + 1))
# 0


We make a small transformation:

$$e^{\frac{2\lambda}{\mu}}e^{\log{\Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu}+1\right)\right)}}$$

$$e^{\frac{2\lambda}{\mu} + \log{\Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu}+1\right)\right)}}$$

Scipy exposes a numerically stable norm.logcdf that we can use:

lambda_ = 1500
mu = 1
x = 1
stats.norm.logcdf(-np.sqrt(lambda_/mu) * (x/mu + 1))
# -3005.26


So now the term in the exponential is much more reasonable:

2 * lambda_ / mu + stats.norm.logcdf(-np.sqrt(lambda_/mu) * (x/mu + 1))
# -5.26

• That seemed like a nice solution but as soon as we have a difference between $x$ and $mu$ greater than $0.3$ we still encounter overflow problems, e.g. lambda_ = 1500 x = 0.70 mu = 1. 2 * lambda_ / mu + norm.logcdf(-np.sqrt(lambda_/mu) * (x/mu + 1))  We obtain $827.3935924745633$ Commented Nov 22, 2020 at 17:09