Repository: scipy/scipy
Pull request: MAINT: fix numerical precision in nct.stats
#### Reference issue
Closes gh-12919

#### What does this implement/fix?
For values of df > 344 the stats method
    of nct would raise runtime warnings and
    become nan unnecessarily.

#### Additional information
There are other known numerical issues with `nct` that ***are not*** fixed by this PR. For example the `cdf` and `pdf` have numerical issues: 
cdf:
https://github.com/scipy/scipy/issues/8344
pdf:
https://github.com/scipy/scipy/issues/7104

It looks like both the `cdf` and `pdf` are (eventually) in fortran so I can't fix. 
Better to get this in and fix a known issue though so opening PR. 

Problem statement:
BUG: numerical precision, use gammaln in nct.mean
possible precision improvement, no bug yet

I was looking at the mean for nct.
I didn't find any problems in my examples with nct mean, but the gamma values were very large 1e60.

using gammaln should improve precision
`np.exp(special.gammaln(df/2.-0.5) - special.gammaln(df/2.))`

```
df = 100
special.gamma(df/2.-0.5), special.gamma(df/2.)
(8.667601843135274e+61, 6.082818640342675e+62)
```

Hints from issue discussion:
ok bug, when it overflows

```
stats.nct.mean(1000, 2)
nan
```

at df=1000, the gamma values are inf

at df=350 it already overflows, mean is nan
AFAIU: c11 should converge to 1   (roughly normal distribution limit of large df)

```
df = 10000
gfac = np.exp(special.gammaln(df/2.-0.5) - special.gammaln(df/2.))
c11 = np.sqrt(df/2.) * gfac
c11
1.000075007811785
``

It looks like the gammaln version is good to precision (6 decimals) at around df = 1e6 or 1e-7. After that c11 bounces around 1.
Hi Josef 

Might be this one:

https://github.com/scipy/scipy/blob/01d8bfb6f239df4ce70c799b9b485b53733c9911/scipy/stats/_continuous_distns.py#L6033

I can see about a test case based on your example here to see if it fixes. 

-Lucas Roberts 

> On Oct 4, 2020, at 11:10 PM, Josef Perktold <notifications@github.com> wrote:
> 
> ﻿
> AFAIU: c11 should converge to 1 (roughly normal distribution limit of large df)
> 
> df = 10000
> gfac = np.exp(special.gammaln(df/2.-0.5) - special.gammaln(df/2.))
> c11 = np.sqrt(df/2.) * gfac
> c11
> 1.000075007811785
> ``
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub, or unsubscribe.

Some numerical confirmation and investigations on my laptop: 

```
>>> from scipy.stats import nct
>>> nct.mean(1000, 2)
/Users/rlucas/scipy-dev/scipy/scipy/stats/_continuous_distns.py:5875: RuntimeWarning: invalid value encountered in true_divide
  gfac = sc.gamma(df/2.-0.5) / sc.gamma(df/2.)
nan
>>> nct.mean(350, 2)
nan
>>> nct.mean(340, 2)
2.004425323005375
>>> nct.mean(345, 2)
nan
>>> nct.mean(343, 2)
2.0043864996470164
>>> nct.mean(344, 2)
0.0
>>> nct.mean(345, 2)
nan
```
So I'm seeing the same thing and on my machine 345 is the limit. 

Now the underlying issue we suspect is this line  

https://github.com/scipy/scipy/blob/01d8bfb6f239df4ce70c799b9b485b53733c9911/scipy/stats/_continuous_distns.py#L6033

so let's confirm that threshold and behavior there

```
>>> from scipy.special import gamma
>>> for df in range(340, 351):
...     gfac = gamma(df/2.-0.5) / gamma(df/2.)
...     print("df", df)
...     print("gfac", gfac)
... 
df 340
gfac 0.07686620227521021
df 341
gfac 0.07675291402619422
df 342
gfac 0.07664012520969488
df 343
gfac 0.07652783216687982
df 344
gfac 0.0
__main__:2: RuntimeWarning: invalid value encountered in double_scalars
df 345
gfac nan
df 346
gfac nan
df 347
gfac nan
df 348
gfac nan
df 349
gfac nan
df 350
gfac nan
``` 
so it looks like this is (a) source of bug. 

EDIT: 

I did a check out to `df=100000` and it looks like the `gammaln` is pretty stable but the convergence to 1 is pretty slow, I'll spare the full output but something like:

```
>>> for df in range(1000, 100001):
...     gfac = exp(gammaln(df/2.-0.5) - gammaln(df/2.))
...     print("gfac", gfac)
...     print("df", df)
...     print(np.sqrt(df/2.) * gfac)
... 
[elided output]
...
gfac 0.004472169496146954
df 100000
1.0000075000285573
```
I'll add in a test for large `df` to `nct` and change the `gfac` line to see if that is all that needs to change.

Commit message: MAINT: fix numerical precision in nct.stats

  * for values of df > 344 the stats method
    of nct would raise runtime warnings and
    become nan unnecessarily.

MAINT: add fix for df < 2 cases to nan

Update scipy/stats/tests/test_distributions.py

Co-authored-by: Warren Weckesser <warren.weckesser@gmail.com>

Update scipy/stats/tests/test_distributions.py

Co-authored-by: Warren Weckesser <warren.weckesser@gmail.com>

bump rtol 1e-13 to 1e-10 after mpmath

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 48891e226fe8..5e4304b8e793 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -5912,12 +5912,12 @@ def _stats(self, df, nc, moments='mv'):
         #
         mu, mu2, g1, g2 = None, None, None, None
 
-        gfac = sc.gamma(df/2.-0.5) / sc.gamma(df/2.)
+        gfac = np.exp(sc.betaln(df/2-0.5, 0.5) - sc.gammaln(0.5))
         c11 = np.sqrt(df/2.) * gfac
-        c20 = df / (df-2.)
+        c20 = np.where(df > 2., df / (df-2.), np.nan)
         c22 = c20 - c11*c11
-        mu = np.where(df > 1, nc*c11, np.inf)
-        mu2 = np.where(df > 2, c22*nc*nc + c20, np.inf)
+        mu = np.where(df > 1, nc*c11, np.nan)
+        mu2 = np.where(df > 2, c22*nc*nc + c20, np.nan)
         if 's' in moments:
             c33t = df * (7.-2.*df) / (df-2.) / (df-3.) + 2.*c11*c11
             c31t = 3.*df / (df-2.) / (df-3.)
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 4cd73e8d943e..6f147b701352 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -3429,14 +3429,39 @@ def test_variance_gh_issue_2401(self):
 
     def test_nct_inf_moments(self):
         # n-th moment of nct only exists for df > n
+        m, v, s, k = stats.nct.stats(df=0.9, nc=0.3, moments='mvsk')
+        assert_equal([m, v, s, k], [np.nan, np.nan, np.nan, np.nan])
+
         m, v, s, k = stats.nct.stats(df=1.9, nc=0.3, moments='mvsk')
         assert_(np.isfinite(m))
-        assert_equal([v, s, k], [np.inf, np.nan, np.nan])
+        assert_equal([v, s, k], [np.nan, np.nan, np.nan])
 
         m, v, s, k = stats.nct.stats(df=3.1, nc=0.3, moments='mvsk')
         assert_(np.isfinite([m, v, s]).all())
         assert_equal(k, np.nan)
 
+    def test_nct_stats_large_df_values(self):
+        # previously gamma function was used which lost precision at df=345
+        # cf. https://github.com/scipy/scipy/issues/12919 for details
+        nct_mean_df_1000 = stats.nct.mean(1000, 2)
+        nct_stats_df_1000 = stats.nct.stats(1000, 2)
+        # These expected values were computed with mpmath. They were also
+        # verified with the Wolfram Alpha expressions:
+        #     Mean[NoncentralStudentTDistribution[1000, 2]]
+        #     Var[NoncentralStudentTDistribution[1000, 2]]
+        expected_stats_df_1000 = [2.0015015641422464, 1.0040115288163005]
+        assert_allclose(nct_mean_df_1000, expected_stats_df_1000[0],
+                        rtol=1e-10)
+        assert_allclose(nct_stats_df_1000, expected_stats_df_1000,
+                        rtol=1e-10)
+        # and a bigger df value
+        nct_mean = stats.nct.mean(100000, 2)
+        nct_stats = stats.nct.stats(100000, 2)
+        # These expected values were computed with mpmath.
+        expected_stats = [2.0000150001562518, 1.0000400011500288]
+        assert_allclose(nct_mean, expected_stats[0], rtol=1e-10)
+        assert_allclose(nct_stats, expected_stats, rtol=1e-10)
+
 
 class TestRice(object):
     def test_rice_zero_b(self):


Code review comments:
diff_hunk: @@ -3429,14 +3429,39 @@ def test_variance_gh_issue_2401(self):
 
     def test_nct_inf_moments(self):
         # n-th moment of nct only exists for df > n
+        m, v, s, k = stats.nct.stats(df=0.9, nc=0.3, moments='mvsk')
+        assert_equal([m, v, s, k], [np.nan, np.nan, np.nan, np.nan])
+
         m, v, s, k = stats.nct.stats(df=1.9, nc=0.3, moments='mvsk')
         assert_(np.isfinite(m))
-        assert_equal([v, s, k], [np.inf, np.nan, np.nan])
+        assert_equal([v, s, k], [np.nan, np.nan, np.nan])
 
         m, v, s, k = stats.nct.stats(df=3.1, nc=0.3, moments='mvsk')
         assert_(np.isfinite([m, v, s]).all())
         assert_equal(k, np.nan)
 
+    def test_nct_stats_large_df_values(self):
+        # previously gamma function was used which lost precision at df=345
+        # cf. https://github.com/scipy/scipy/issues/12919 for details
+        nct_mean_df_1000 = stats.nct.mean(1000, 2)
+        nct_stats_df_1000 = stats.nct.stats(1000, 2)
+        # These expected values were computed with mpmath. They were also
+        # verified with the Wolfram Alpha expressions:
+        #     Mean[NoncentralStudentTDistribution[1000, 2]]
+        #     Var[NoncentralStudentTDistribution[1000, 2]]
+        expected_stats_df_1000 = [2.0015015641422464, 1.0040115288163005]
+        assert_allclose(nct_mean_df_1000, expected_stats_df_1000[0],
+                        rtol=1e-10)
+        assert_allclose(nct_stats_df_1000, expected_stats_df_1000,
+                        rtol=1e-10)
+        # and a bigger df value
+        nct_mean = stats.nct.mean(100000, 2)
+        nct_stats = stats.nct.stats(100000, 2)
+        # These expected values were computed with mpmath.
+        expected_stats = [2.0000150001562518, 1.0000400011500288]
+        assert_allclose(nct_mean, expected_stats[0], rtol=1e-10)
+        assert_allclose(nct_stats, expected_stats, rtol=1e-10)
path: scipy/stats/tests/test_distributions.py
text: @user1:
This is the test that is failing, which means it is the second component of `nct_stats` (i.e. the variance) that has the larger relative error.  The tests should pass if we bump up `rtol` to 1e-9.
```suggestion
        assert_allclose(nct_stats, expected_stats, rtol=1e-9)
```

@author:
> This is the test that is failing, which means it is the second component of `nct_stats` (i.e. the variance) that has the larger relative error. The tests should pass if we bump up `rtol` to 1e-9.

Ah, ok. This explains then, the previous assertion failure was the line just prior that tests on the mean. 
This means a bump to `rtol=1e-10` fixed that assertion. 

I committed the suggestion on the `rtol=1e-9` here for the variance, hopefully that is the last of the precision failing tests.