ENH: special: use boost in nctdtr

[dschmitz89]

Reference issue

Closes #8344

What does this implement/fix?

Uses boost for nctdtr. The PR is modeled after #21505. The trickiest part was the mpmath reference implementation,.

Additional information

I did not migrate the noncentral t quantile nctdtrit yet as all default tests (no xslow etc.) ran successfully locally and I have to adjust the size of my PRs to the sleeping cycles of my little one, so trying to keep the small.

[dschmitz89]

This is the test of the parameter range in the linked issue.

[fancidev]

Looks good to me! Just some small comments below.

[dschmitz89]

@josef-pkt @bashtage : would migrating nctdtr to boost potentially cause issues for statsmodels? Like with ncfdtr, for the vast majority of cases we get much better precision.

In case you have roundtrip tests for nctdtrit, nctdtridf or nctdtrinc, they might fail for extreme inputs though. Don't worry though, scipy.special contains its own low level root finder now, so whatever CDF we use, we can invert it with respect to any parameter we want in principle.

[josef-pkt]

I switched away from scipy stats distributions for both nct and ncf to use scipy special
statsmodels/statsmodels#8659

AFAIR, AFAICS for nct it was only a runtime warning that showed up during rootfinding in sample size computation.
It did not affect the final results. So, I don't know which corner/extreme case evaluated during the rootfinding caused the Runtime Warning.

(Because I always use older versions of scipy, lagging 1 to 3 years behind, I cannot investigate new problems when they show up.)

[dschmitz89]

I switched away from scipy stats distributions for both nct and ncf to use scipy special statsmodels/statsmodels#8659

AFAIR, AFAICS for nct it was only a runtime warning that showed up during rootfinding in sample size computation. It did not affect the final results. So, I don't know which corner/extreme case evaluated during the rootfinding caused the Runtime Warning.

(Because I always use older versions of scipy, lagging 1 to 3 years behind, I cannot investigate new problems when they show up.)

Looking at statsmodels/statsmodels#8624, it would be interesting to know which parameters caused Boost's noncentral t CDF to go wild.

[josef-pkt]

I have now scipy 0.10

>>> stats.nct.cdf(4.30265273, 2., 0.978343)
Traceback (most recent call last):
  File "<pyshell#15>", line 1, in <module>
    stats.nct.cdf(4.30265273, 2., 0.978343)
  File "...\WPy64-31140\python-3.11.4.amd64\Lib\site-packages\scipy\stats\_distn_infrastructure.py", line 2186, in cdf
    place(output, cond, self._cdf(*goodargs))
  File "...\WPy64-31140\python-3.11.4.amd64\Lib\site-packages\scipy\stats\_continuous_distns.py", line 6826, in _cdf
    return np.clip(_boost._nct_cdf(x, df, nc), 0, 1)
RuntimeWarning: invalid value encountered in _nct_cdf

this is just looking at the values in the debugger of reported case in statsmodels issue
(I change the statsmodels code to use scipy.stats nct.cdf)

update
If I don't force the warning to error, then I get the number after the warning

>>> stats.nct.cdf(4.30265273, 2., 0.978343)

Warning (from warnings module):
  File "...\WPy64-31140\python-3.11.4.amd64\Lib\site-packages\scipy\stats\_continuous_distns.py", line 6826
    return np.clip(_boost._nct_cdf(x, df, nc), 0, 1)
RuntimeWarning: invalid value encountered in _nct_cdf
0.9107619358152059

correction: initial warning in statsmodels example is in sf, but that raises the same warning as cdf
stats.nct.sf(4.30265273, 2., 0.978343)

...
    return np.clip(_boost._nct_sf(x, df, nc), 0, 1)
RuntimeWarning: invalid value encountered in _nct_sf

[bashtage]

I don't recall the specifics TBH. It was failing in pre-release testing which led to the previous bug report.

[dschmitz89]

stats.nct.cdf(4.30265273, 2., 0.978343)

Both warnings do not appear anymore with a recent scipy version. Maybe you could just suppress it.

This PR is ready from my side.

[josef-pkt]

Thanks,
We can keep using the scipy special version of nct in statsmodels for now, then users that upgrade to the next scipy will get the fixed nct boost version.

[dschmitz89]

@steppi: if you could provide an example of how to use the new rootfinder with boost functions, I could tackle migrating a lot of cdflib functions over to boost over the next months piece by piece. Such PRs are not really popular as they are tedious but they fit my current time budget for individual PRs reasonably well.

[steppi]

Why are we removing these test cases? Were they failing because the old cdflib reference is inaccurate for some of them? The reference values could be replaced with mpmath then. This seems like a good opportunity to remove the xfail_on_32bit and show that they now pass.

[dschmitz89]

Yes, the old reference values were very inaccurate for some cases. Allright, I will integrate them in the new tests against the mpmath reference. I might make a subjective selection though if that's okay with you.

[dschmitz89]

For high noncentrality values of 38 such as in the test originally, the mpmath reference computes CDF values of the order of 1e-310. There also boost fails to get close or returns nan, so I dropped these for noncentrality of 15 instead.

[steppi]

Cases like this where the reference output is subnormal are why it makes sense to use an atol in tests. In cases like this, it's fine to return zero and signal underflow through sf_error. Is Boost returning nan where the previous method used to return 0? I'd count that as a bug.

[bashtage]

Fully agree that the nan feels like a bug.

[dschmitz89]

@steppi : I do , am just short on time at the moment. Not sure I can make it till the 1.15 release as I am also going on holiday at the end of the month. I want to spend my scipy activities until then mostly on review. If you can/want to take over the last todos here, that would be much appreciated as I think this PR is almost ready to merge.

Edit: I can also file an issue upstream about the nan values but I doubt that they will be fixed in a short time frame.

[dschmitz89]

@steppi : I just had another look at the test cases for nc=38. Even for small values within the normal floating point range, we get high discrepancies.

For df, nc, t = 0.98, 38, 1.5 the mpmath reference computes 1.546906350750292e-295 but boost results in 2.5919953604832226e-97. Would you prefer to get the test to pass with something like atol=1e-90? I would simply ignore these edge cases for now to be honest. In the parameter ranges that matter most we get much better results with boost.

[steppi]

It looks like you set mp.dps too low to evaluate the sum. I set mp.dps=500 and got the result 2.591995360483094e-97 from your reference implementation. Much closer to Boost. If you're evaluating a sum with mpmath, you may need to set dps pretty high.

[dschmitz89]

Ouch, thanks for spotting this. With mp.dps=500, the reference calculation did not finish within 10 minutes on my laptop though. It has performance issues at the moment but this is crazy. How long does it take on your machine? Might try out my Mac instead.

[dschmitz89]

With 400 significant digits I also get 2.591995360483094e-97 within 12 minutes of computation time.

[steppi]

@steppi: if you could provide an example of how to use the new rootfinder with boost functions, I could tackle migrating a lot of cdflib functions over to boost over the next months piece by piece. Such PRs are not really popular as they are tedious but they fit my current time budget for individual PRs reasonably well.

Sure, that would be great. I'll do one when I get bandwidth for it, and will ping you on the PR. This would likely be after the xsf migration, which I hope to complete before we branch.

[dschmitz89]

Allright, stars aligned and I had a bit of time yesterday.

The nan value originates from nctdtr(980, 38, 1.5), here boost computes a tiny negative value of app. -1e-70. In the wrapper, this is set to nan. I worked around it in the test by ignoring it. We could theoretically also clip for tiny negative values to zero but who knows if that is the better option for other cases.

There is now one failing parameter combination for 32bit CPUs. Let's skip the test again on 32bit then?

@steppi : let me know if you have a strong opinion here about any of the above points, otherwise I will add the test skip.

[steppi]

@dschmitz89, instead of removing the failing test case, mark it as xfail with a reason, and then do a conditional xfail for only the one failing case on 32bit.

Syntax is, instead of having a test case

[0.98, -3.8, 0.15, 0.9999528361700505, 1e-13], in the pytest.markparametrize` list, put

            pytest.param(
                [0.98, -3.8, 0.15, 0.9999528361700505, 1e-13],
                marks=pytest.mark.xfail(reason="Bug in underlying Boost math implementation"),
            )

for the failing test case on 32bit, your pytest.param will look like

pytest.param(
    test_case,  # Subsitute in your actual test case
    marks=pytest.mark.xfail(condition=sys.maxsize < 2**32, reason="Fails on 32 bit.")
)

I don't think it's a good idea to just get rid of failing test cases, or skip the whole test on 32bit because of only one failing case.

[dschmitz89]

Interesting that the 32bit job succeeds now. Maybe it had an issue in the way I handled the nan case which is now xfailed. Will push the latest stylistic updates without rerunning CI.

[dschmitz89]

The green CI results do not show up here anymore due to my last commit: for the reviewer, they can easily be found here.