BUG: stats.nct.pdf inconsistent behavior when compared to MATLAB and R

[Martin2Constant]

Describe your issue.

Hello,

When computing stats.nctpdf with large t and nc values, the function returns 'inf'.
This behavior is inconsistent with both MATLAB and R. This "bug" appeared with version scipy 1.10 (behavior was consistent before).
I'm not sure whether it's truly a bug or whether MATLAB and R are both wrong in not returning inf.

The following R code (equivalent to the Python code) returns 0.0043403 instead of 'inf' (R version 4.3.1):

t <- 8.45316107575123
N <- 25
df <- N-1
safe_int <- .9999
safe_range <- t / sqrt(N) + c(-1, 1) * qt(1 - (1 - safe_int) / 2, df) / sqrt(N)

dt(t, df, ncp=safe_range[[2]]*sqrt(N))

Same value is returned in MATLAB (R2023b):

t = 8.45316107575123;
N = 25;
df = N-1;
safe_int = .9999;
safe_range = t / sqrt(N) + [-1, 1] * tinv(1 - (1 - safe_int) / 2, df) / sqrt(N);

nctpdf(t, df, safe_range(2)*sqrt(N))

Reproducing Code Example

import numpy as np
from scipy import stats

t = 8.45316107575123
N = 25
df = N-1
safe_int = .9999
safe_range = t / np.sqrt(N) + np.array([-1, 1]) * stats.t.ppf(1 - (1 - safe_int) / 2, df=df) / np.sqrt(N)

# SciPy 1.9.3 returns 0.0043403; SciPy 1.10 and above return Inf
stats.nct.pdf(t, df, nc=safe_range[1]*np.sqrt(N))

Error message

N/A

SciPy/NumPy/Python version and system information

1.11.3 1.26.0 sys.version_info(major=3, minor=10, micro=12, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: C:/Users/constama/miniforge3/Library/include
    lib directory: C:/Users/constama/miniforge3/Library/lib
    name: blas
    openblas configuration: unknown
    pc file directory: C:\bld\scipy-split_1696467770591\_h_env\Library\lib\pkgconfig
    version: 3.9.0
  lapack:
    detection method: pkgconfig
    found: true
    include directory: C:/Users/constama/miniforge3/Library/include
    lib directory: C:/Users/constama/miniforge3/Library/lib
    name: lapack
    openblas configuration: unknown
    pc file directory: C:\bld\scipy-split_1696467770591\_h_env\Library\lib\pkgconfig
    version: 3.9.0
  pybind11:
    detection method: pkgconfig
    include directory: C:/Users/constama/miniforge3/Library/include
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    commands: clang-cl
    linker: lld-link
    name: clang-cl
    version: 17.0.0
  c++:
    commands: clang-cl
    linker: lld-link
    name: clang-cl
    version: 17.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.36
  fortran:
    commands: flang-new
    linker: lld-link
    name: flang
    version: 17.0.0
  pythran:
    include directory: ..\..\_h_env\lib\site-packages\pythran
    version: 0.14.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
Python Information:
  path: C:\bld\scipy-split_1696467770591\_h_env\python.exe
  version: '3.10'

[WarrenWeckesser]

I'm not sure whether it's truly a bug or whether MATLAB and R are both wrong in not returning inf.

Thanks for reporting the issue. It is a bug, and it is still present in the development version of SciPy. Here's a simplified example, with integer parameter values so it is easy to compare to the result computed by Wolfram Alpha:

In [28]: import scipy

In [29]: scipy.__version__
Out[29]: '1.12.0.dev0+1813.0a3dce2'

In [30]: from scipy.stats import nct

In [31]: nct.pdf(8, 24, 13)
Out[31]: inf

With Wolfram Alpha, PDF[NoncentralStudentTDistribution[24, 13], 8] gives 0.0017644351379120543496495442757378120219461588078650315935117814288005136461742269261929088489060060675999625476400913578446784114357283664002733506219393887469961463578572949425885674236825307738127452 (after clicking on "More digits" a couple times).

[josef-pkt]

too much boost ?

>>> from scipy import stats
>>> stats.nct.pdf(8, 24, 13)
0.0017644351379121064
>>> stats.nct.pdf(8, 24, 13) - .0017644351379120543496
5.204170427930421e-17

>>> import scipy
>>> scipy.__version__
'1.7.3'

[mdhaber]

Boost isn't used by nct.pdf currently.
https://github.com/scipy/scipy/blob/main/scipy/stats/_continuous_distns.py#L7631-L7709

[WarrenWeckesser]

The development version of boost/math gives the correct result in this case. noncentralt_pdf.cpp takes x, df and nc from the command line and evaluates the PDF using boost's non_central_t distribution:

$ ./noncentralt_pdf 8 24 13
x  = 8
df = 24
nc = 13
p  =  1.7644351379120544e-03

As noted in the comment in the SciPy code, the boost implementation has issues in the left tail.

[mdhaber]

@WarrenWeckesser are you willing to open an issue with Boost about the left tail of the pdf? It would be nice if we could switch to Boost's implementation after all rather than stitching ours and theirs together. Here's a demonstration in SciPy (using Boost for the PDF):

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
dist = stats.nct(8, 8.5)
x = np.linspace(-3, 5)
plt.semilogy(x, dist.pdf(x))
plt.show()

Or maybe in the meantime just going with Boost is the lesser of two evils.

[josef-pkt]

one possible candidate for change compared to my version is
#17302 boost

(because pdf code itself is unchanged0

[WarrenWeckesser]

Thanks @josef-pkt, that is in fact the source of the problem reported here. When the b parameter of hyp1f1 is 1.5 and the third parameter (x) is greater than 50, hyp1f1 returns inf:

In [17]: import scipy

In [18]: scipy.__version__
Out[18]: '1.12.0.dev0+1814.166e1f2'

In [19]: from scipy.special import hyp1f1

In [20]: hyp1f1(13.0, [1.499999, 1.5, 1.500001], 61.0).tolist()
Out[20]: [1.3550915047757008e+39, inf, 1.3550800368458125e+39]

The problem is that the underlying boost function hypergeometric_1F1 throws an exception when it shouldn't. I reported the boost issue here: boostorg/math#1034.

I'll see if we can work around the boost bug in our wrapper of hypergeometric_1F1. This needs to be fixed regardless of whether or not we replace the current implementation of the PDF with the boost noncentral t PDF function.

The problem with the PDF calculation in the left tail is a separate issue. @mdhaber, if you already have the evidence, go ahead and file a boost issue. I probably won't get to that in the immediate future.

[mdhaber]

I ask because I don't ever work with C++. We'll see if they'll take the report in Python form.

[jzmaddock]

There is a fix for the left tail issue in the works at the Boost end now, however, the accuracy is still poor and I see no easy fix at present. Would I be correct in thinking that a handful of digits correct is "good enough" for most stats usages?

[mdhaber]

Probably (although we are regularly surprised at what is requested). In any case, I think that fixing that discontinuity would be enough of an improvement to switch to the Boost implementation so we can fix the garbage output reported in this issue.

Update: when we do fix this, check the cases in gh-19450.

[boricua-hg]

Can I ad that if one wants to use scipy to develop methods used in validation then accuracy is important? Was a user of Minitab in previous life in quality and developing sampling plans was part of the job. Accuracy similar to Minitab would be preferable.

[boersmamarcel]

@jzmaddock is this now already part of the SciPy release? I am running into the same problems. Running SciPy 1.12.0.

[jzmaddock]

I can't comment on SciPy, but this would have been in the last Boost release (1.84), so it probably depends what Boost release is installed when SciPy is built?

[lucascolley]

SciPy was recently updated to use 1.83 on main, so would need an update before the next release. It sounds like that isn't trivial and introduces some errors which would need resolving.

[jzmaddock]

Anything on our side that we need to look at?

[abdrysdale]

From looking at the source code, I managed to fix this myself by converting the calculation into log space before converting back. The problem with the scipy implementation (at least for me) was that the numbers got exceptionally large before being reduced down.

I've written my own implementation for a project that I'm using Pytorch but the idea is the same

Pi = torch.acos(torch.zeros(1)) * 2

def nct_pdf(x, df, nc, norm=True):
    """Non-central t distribution probability density function.

    There are problems with the scipy implementation so this function fixes
    those by performing most of the calculations in log space.
    Most of the function is also implemented in pytorch.
    Maths taken from: https://en.wikipedia.org/wiki/Noncentral_t-distribution#Probability_density_function

    Args:
        x (torch tensor) : Values to evaluate the pdf at.
        df (int) : Degrees of freedom.
        nc (float) : Non-centrallity parameter.

    Returns:
        pdf (torch.tensor) : Probability density at each x.
    """
    n = torch.tensor(df)
    nc = torch.tensor(nc)

    x2 = x ** 2
    mu2 = nc ** 2
    mu2x2 = mu2 * x2
    _2 = torch.tensor(2)

    # Student T(x; mu = 0)
    lgamma2 = torch.lgamma(n / 2)
    lgamma12 = torch.lgamma((n + 1) / 2)
    stmu_z = (
        lgamma12 - lgamma2 - ((n + 1)/2) * torch.log(1 + x2/n)
        - 0.5 * torch.log(n * Pi) - mu2/2
    )

    # A_v(x; mu) and B_v(x; mu) terms
    zterm = mu2x2 / (2 * (x2 + n))
    Av = loghyp1f1((df + 1)/2, .5, zterm)
    Bv = (
        torch.log(torch.sqrt(_2) * nc * x / torch.sqrt(x2 + n))
        + torch.lgamma(n/2 + 1) - lgamma12
        + loghyp1f1(df/2 + 1, 1.5, zterm)
    )

    pdf = torch.exp(stmu_z + Av) + torch.exp(stmu_z + Bv)
    pdf = torch.nan_to_num(pdf)
    if norm:
        pdf /= torch.sum(pdf)
        pdf = torch.nan_to_num(pdf)
    return torch.clamp(pdf, 0, 1).to(x.dtype)

where the function loghyp1f1 is the logarithm of the confluent hypergeometric function that I've written in Fortran and that I'm calling from Python.

Happy to submit a pull request if this would be useful.

[lucascolley]

I think that I get 2 extra errors with boost 1.84. Not sure which end they are due to.

scipy/stats/tests/test_continuous_basic.py:170: in test_cont_basic
    check_sf_isf(distfn, arg, distname)
        arg        = ()
        distfn     = <scipy.stats._continuous_distns.wald_gen object at 0x119987a00>
        distname   = 'wald'
        m          = 1.0
        n_fit_samples = 200
        rng        = RandomState(MT19937) at 0x11E5B2F40
        rvs        = array([2.00676157, 0.29004979, 0.62849545, 1.6186909 , 0.3169379 ,
       0.24333655, 1.62317229, 0.62463867, 3.540782...25, 0.33951168, 2.29954848, 0.79173576, 0.32611827,
       2.27908213, 0.92946228, 1.14471458, 0.27071012, 1.77247785])
        sn         = 500
        v          = 1.0
scipy/stats/tests/test_continuous_basic.py:597: in check_sf_isf
    npt.assert_almost_equal(distfn.sf(distfn.isf([0.1, 0.5, 0.9], *arg), *arg),
        arg        = ()
        distfn     = <scipy.stats._continuous_distns.wald_gen object at 0x119987a00>
        msg        = 'wald'
scipy/stats/_distn_infrastructure.py:2164: in sf
    place(output, cond, self._sf(*goodargs))
        _a         = 0.0
        _b         = inf
        args       = ()
        cond       = array([1, 1, 1])
        cond0      = 1
        cond1      = array([ True,  True,  True])
        cond2      = array([0, 0, 0])
        dtyp       = dtype('float64')
        goodargs   = [array([4.47482874e+248, 6.75841306e-001, 2.37624709e-001])]
        kwds       = {}
        loc        = array(0)
        output     = array([0., 0., 0.])
        scale      = array(1)
        self       = <scipy.stats._continuous_distns.wald_gen object at 0x119987a00>
        x          = array([4.47482874e+248, 6.75841306e-001, 2.37624709e-001])
scipy/stats/_continuous_distns.py:10743: in _sf
    return invgauss._sf(x, 1.0)
        self       = <scipy.stats._continuous_distns.wald_gen object at 0x119987a00>
        x          = array([4.47482874e+248, 6.75841306e-001, 2.37624709e-001])
scipy/stats/_continuous_distns.py:4784: in _sf
    return np.exp(self._logsf(x, mu))
        mu         = 1.0
        self       = <scipy.stats._continuous_distns.invgauss_gen object at 0x118fdf7c0>
        x          = array([4.47482874e+248, 6.75841306e-001, 2.37624709e-001])
scipy/stats/_continuous_distns.py:4781: in _logsf
    return a + np.log1p(-np.exp(b - a))
E   RuntimeWarning: divide by zero encountered in log1p
        a          = array([-2.23741437e+248, -4.25683513e-001, -6.07213810e-002])
        b          = array([-2.23741437e+248, -1.87520797e+000, -3.19210227e+000])
        fac        = array([4.72728505e-125, 1.21640343e+000, 2.05141819e+000])
        mu         = 1.0
        self       = <scipy.stats._continuous_distns.invgauss_gen object at 0x118fdf7c0>
        x          = array([4.47482874e+248, 6.75841306e-001, 2.37624709e-001])
---
scipy/stats/tests/test_fast_gen_inversion.py:118: in test_rvs_and_ppf
    assert_allclose(rng1.ppf(q), rng2.ppf(q), atol=1e-10)
        args       = ()
        distname   = 'wald'
        q          = [0.001, 0.1, 0.5, 0.9, 0.999]
        rng1       = <scipy.stats._distn_infrastructure.rv_continuous_frozen object at 0x12000d7e0>
        rng2       = <scipy.stats._sampling.FastGeneratorInversion object at 0x12000dd50>
        rvs1       = array([0.49415957, 0.08335714, 0.32718863, 0.68736644, 1.00888463,
       0.21064594, 0.27862896, 3.03444895, 0.356945...43, 1.01103814, 0.22827809, 0.4331922 , 0.31824336,
       2.24261899, 0.24624813, 0.63768449, 0.66763203, 1.55694817])
        rvs2       = array([1.13945461, 2.10754385, 0.44018833, 1.04186099, 0.15792175,
       1.18501888, 0.60517616, 0.12677009, 0.557995...38, 0.16108085, 0.2168502 , 3.36257858, 0.33278517,
       2.67267979, 2.19542446, 2.64032752, 0.28419598, 0.75260723])
        urng       = Generator(PCG64) at 0x12085C9E0
../../../../../../mambaforge/envs/scipy-dev/lib/python3.10/contextlib.py:79: in inner
    return func(*args, **kwds)
E   AssertionError: 
E   Not equal to tolerance rtol=1e-07, atol=1e-10
E   
E   Mismatched elements: 1 / 5 (20%)
E   Max absolute difference: 4.47482874e+248
E   Max relative difference: 2.08808116e+248
E    x: array([7.921848e-002, 2.376247e-001, 6.758413e-001, 4.474829e+248,
E          8.354865e+000])
E    y: array([0.079218, 0.237625, 0.675841, 2.143034, 8.354865])
        args       = (<function assert_allclose.<locals>.compare at 0x120206950>, array([7.92184778e-002, 2.37624709e-001, 6.75841306e-001, 4.47482874e+248,
       8.35486493e+000]), array([0.07921848, 0.23762471, 0.67584131, 2.14303391, 8.35486495]))
        func       = <function assert_array_compare at 0x116558e50>
        kwds       = {'equal_nan': True, 'err_msg': '', 'header': 'Not equal to tolerance rtol=1e-07, atol=1e-10', 'verbose': True}
        self       = <contextlib._GeneratorContextManager object at 0x11652f130>

[jzmaddock]

From looking at the source code, I managed to fix this myself by converting the calculation into log space before converting back. The problem with the scipy implementation (at least for me) was that the numbers got exceptionally large before being reduced down.

If SciPy is just calling the Boost implementation then we shouldn't spuriously overflow - if we do it's a bug and we'll fix that.

I think that I get 2 extra errors with boost 1.84. Not sure which end they are due to.

I fear someone will have to explain that output to this non-Python person ;)

[abdrysdale]

From looking at the source code, I managed to fix this myself by converting the calculation into log space before converting back. The problem with the scipy implementation (at least for me) was that the numbers got exceptionally large before being reduced down.

If SciPy is just calling the Boost implementation then we shouldn't spuriously overflow - if we do it's a bug and we'll fix that.

Is scipy just calling boost? I've found some nct pdf code here - if this is being used it might benefit from performing the bulk of the calculations in log space.

[dschmitz89]

From looking at the source code, I managed to fix this myself by converting the calculation into log space before converting back. The problem with the scipy implementation (at least for me) was that the numbers got exceptionally large before being reduced down.

If SciPy is just calling the Boost implementation then we shouldn't spuriously overflow - if we do it's a bug and we'll fix that.

I think that I get 2 extra errors with boost 1.84. Not sure which end they are due to.

I fear someone will have to explain that output to this non-Python person ;)

The new errors with BooostMath 1.84 are unrelated to this issue. They originate from Boost's implementation of the inverse survival function of the inverse gaussian distribution, see here. Wald is a special case of the inverse gaussian distribution for $\mu=1$. The reported errors look like the inverse survival function blows up for $q=0.1$ or $q=0.9$ on some Mac hardware.

Side note: the RunTimeError in logsf does not look good either, will see if something can be done there.

[jzmaddock]

The new errors with BooostMath 1.84 are unrelated to this issue. They originate from Boost's implementation of the inverse survival function of the inverse gaussian distribution, see here. Wald is a special case of the inverse gaussian distribution for . The reported errors look like the inverse survival function blows up for or on some Mac hardware.

Hmmm, I don't have a mac here, but I have tried to reproduce and completely failed on Windows/MSVC, @mborland are you able to try on your Mac?

Here's the test program I used, it tests every representable q value from 0.09999999 to 0.10000001, and takes a good couple of hours to run even in release mode. There is a bit of "wobble" in the quantiles returned (up to 9ULP on my machine), but that's to be expected given that (a) the root finder stops once it's "close enough" and (b) there may be multiple abscissa values that are equally correct. Which is to say, the quantile doesn't quite always monotonically decrease for increasing q.

#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false

#include <iostream>
#include <boost/math/distributions.hpp>

void print_error(double val, int distance)
{
   std::cout << "New high error of " << distance << " ULP at " << std::setprecision(std::numeric_limits<double>::max_digits10) << val << std::endl;
}

int main()
{
   boost::math::inverse_gaussian g;

   double q = 0.09999999;

   double p1_last = 5;
   double p2_last = 5;

   int distance_1 = 0;
   int distance_2 = 0;

   std::cout << std::setprecision(std::numeric_limits<double>::max_digits10) <<
      quantile(complement(g, 0.09)) << " " << quantile(complement(g, 0.1)) << " " <<
      quantile(complement(g, 0.11)) << std::endl;

   while (q < 0.10000001)
   {
      q = boost::math::float_next(q);

      double p1 = quantile(g, 1 - q);
      double p2 = quantile(complement(g, q));

      if (!boost::math::isnormal(p1))
         abort();
      if (!boost::math::isnormal(p2))
         abort();

      int d = boost::math::float_distance(p1_last, p1);
      if (d > distance_1)
      {
         print_error(q, d);
         distance_1 = d;
      }
      d = boost::math::float_distance(p2_last, p2);
      if (d > distance_2)
      {
         print_error(q, d);
         distance_2 = d;
      }
      p1_last = p1;
      p2_last = p2;
   }

   return 0;
}

[jzmaddock]

Update: that was tested against develop, there was an old root finding bug which resurfaced in 1.84 causing excessive iterations to be used, I need to double check that...

[jzmaddock]

This is our bad: it's a resurfacing of boostorg/math#184 as a result of accepting a PR we shouldn't have :( The issue is now better tested and fixed again in 1.85 which will be out shortly (I hope). Apologies all round.

[dschmitz89]

This is our bad: it's a resurfacing of boostorg/math#184 as a result of accepting a PR we shouldn't have :( The issue is now better tested and fixed again in 1.85 which will be out shortly (I hope). Apologies all round.

Thanks for taking a look so quickly!

[lucascolley]

thanks a lot @jzmaddock ! I think branching for the next SciPy version is happening very soon, so 1.85 may not make it in this time around - it should make it for the next minor release though, which by the sounds of it should help address this issue. Hope that answers your question @boersmamarcel !

[mdhaber]

This also seems to be fixed in main, so closing.

from scipy import stats
from scipy.special import hyp1f1

print(stats.nct.pdf(8, 24, 13))  # 0.0017644351379121057
print(hyp1f1(13.0, [1.499999, 1.5, 1.500001], 61.0))  # [1.35509151e+39 1.35508577e+39 1.35508004e+39]

Thanks to everyone who worked on this!