scipy: `stats.norminvgauss` incorrect implementation?

Hi,

Mentioned in the latest release of TensorFlow Probability, it has been noted that scipy.stats.norminvgauss has not been implemented in the same way as R or the original mathematical paper.

It seems the loc and scale parameters should actually be mu and delta. Is there a way to update this so that non-standard cases (other than loc=0, scale =1) make mathematical sense?

Note that `scipy.stats.norminvgauss` implements the distribution as a
      location-scale family. However, in the original paper, and other
      implementations (such as R) do not implement it this way. Thus the
      `scale` parameters here and scipy don't match unless `scale = 1`.
  Warning: As mentioned above, this distribution is __not__ a location-scale
  family.

Link here: https://github.com/tensorflow/probability/blob/main/tensorflow_probability/python/distributions/normal_inverse_gaussian.py

Docs: https://www.tensorflow.org/probability/api_docs/python/tfp/distributions/NormalInverseGaussian

About this issue

Original URL
State: closed
Created 3 years ago
Comments: 16 (9 by maintainers)

Most upvoted comments

@mdhaber Thanks! That works really well!

Nishan-Pradhan on Aug 31, 2021

Very close, but don’t use norminvg to calculate the PPF and ISF. The functions will only be evaluated at x = tol/10, so you can code in correct values yourself. It is not very sensitive to these values being super accurate - e.g., for PPF, return any value x at which norminvg.cdf(x) is very small (ideally tol/10).

Something like

class CustomNIG:

    def pdf(self, x):
        return norminvg.pdf(x)

    def cdf(self, x):
        return norminvg.cdf(x)

    def ppf(self, x):
        return -2  # norminvg.cdf(-2) ~ 0

    def isf(self, x):
        return 2  # norminvg.sf(2) ~ 0


dist = CustomNIG()
norminvg_fast = stats.NumericalInverseHermite(dist, tol=1e-6)

ppf = norminvg_fast.ppf(data)

This will make NumericalInverseHermite form an approximation H of the inverse of norminvg.cdf that is valid between -2 and 2. The round-trip tolerance (i.e. abs(x - norminvg.cdf(H(x)))) will be controlled at mesh points and midpoints between mesh points to tol=1e-6.

mdhaber on Aug 26, 2021

I can confirm the slow speed in the development version of SciPy:

import scipy.stats as stats
import numpy as np
import timeit

data  = np.random.rand(100,1)
norminvg = stats.norminvgauss(a=1.091, b=-0.649, loc=0.124, scale=0.168)
ppf  = norminvg.ppf(data)

%timeit -n 1 norminvg.ppf(data)
1.4 s ± 5.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit -n 1 stats.student.ppf(data, 10)
Traceback (most recent call last):

%timeit -n 1 stats.t.ppf(data, 10)
339 µs ± 54.3 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit -n 1 stats.norm.ppf(data)
94 µs ± 39.2 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit -n 1 stats.invgamma.ppf(data, 5)
141 µs ± 40.7 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

But please consider upgrading; we can’t support such an old version.

If you upgrade to 1.7, you may find NumericalInverseHermite useful:

>>> norminvg_fast = stats.NumericalInverseHermite(norminvg, tol=1e-8)  # this takes ~1.5s, but...
>>> %timeit -n 1 norminvg_fast.ppf(data)
The slowest run took 4.65 times longer than the fastest. This could mean that an intermediate result is being cached.
35.7 µs ± 26.5 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

>>> np.max(np.abs(norminvg.ppf(data) - norminvg_fast.ppf(data)))
6.982312505954269e-09

mdhaber on Aug 13, 2021

In TF, if you use loc = mu, tailweight = alpha, skewness = beta, scale = sigma to evaluate the CDF, I think this corresponds to stats.norminvgauss(alpha*sigma, beta*sigma, loc=mu, scale=sigma).cdf(x) in SciPy.

chrisb83 on Aug 10, 2021