2
$\begingroup$

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"])

enter image description here

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"])

enter image description here

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
    numerator = numerator_f1*(numerator_f2_add1_f1*numerator_f2_add1_f2*numerator_f2_add1_f3 + numerator_f2_add2)
    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,

$\endgroup$
3
  • $\begingroup$ (x/mu -1)) * np.exp(2*lamb/mu should be a +, no? $\endgroup$ – Cam.Davidson.Pilon Nov 22 '20 at 13:36
  • $\begingroup$ what's the order of magnitude of $\mu$? $\endgroup$ – Cam.Davidson.Pilon Nov 22 '20 at 13:39
  • $\begingroup$ @Cam.Davidson.Pilon You're totally right, I have no idea how I missed that. $\mu$ is always $0.6<\mu<1.2$ $\endgroup$ – ИванКарамазов Nov 22 '20 at 13:51
1
$\begingroup$
  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
$\endgroup$
1
  • 1
    $\begingroup$ 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$ $\endgroup$ – ИванКарамазов Nov 22 '20 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.