Repository: scipy/scipy
Pull request: BUG: stats: handle infinite `df` in `t` distribution
#### Reference issue

fixes #14777

#### What does this implement/fix?

`t` distribution returned garbage values when `df` was infinity. It is known
that the `t` distribution limits to the normal distribution when `df` approaches
infinity. It is much better to handle infinite `df` by branching to the normal
distribution instead of returning wrong values.

appropriate tests have been added for the fix.

#### Additional information

Problem statement:
BUG: Wrong limit and no warning in stats.t for df=np.inf
### Describe your issue.

Scipy's stats module for Student's t distribution does not have the correct limit of a Gaussian distribution if the degree-of-freedom (`df`) parameter is set to `np.inf`. Most importantly, Scipy does not print a warning nor does it raise an error, instead it simply yields a completely wrong result.

While I think it would make sense to branch for `df=np.inf`, it certainly does not make sense to let the user plug in values that produce bogus results.

### Reproducing Code Example

```python
In [1]: import numpy as np

In [2]: from scipy.stats import norm as stats_norm

In [3]: from scipy.stats import t as stats_t

In [4]: stats_t.sf(loc=3., df=99999., x=2., scale=1.)
Out[4]: 0.8413435362058463

In [5]: stats_t.sf(loc=3., df=np.inf, x=2., scale=1.)
Out[5]: 0.5
```


### Error message

```shell
None
```


### SciPy/NumPy/Python version information

1.7.1 1.21.2 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
Thanks for reporting, @Edenhofer!

> ... I think it would make sense to branch for `df=np.inf` ...

I agree. In fact, that's exactly what `multivariate_t` does right now:

```py
In [1]: from scipy.stats import norm, multivariate_t                                                                                                                                                       

In [2]: norm.pdf(x=2, loc=3, scale=1)                                                                                                                                                                      
Out[2]: 0.24197072451914337

In [3]: multivariate_t.pdf(x=2, df=float('inf'), loc=3)                                                                                                                                                    
Out[3]: 0.24197072451914337
```

Commit message: MAINT: stats: use lazy functions only

Patch to review:
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 4d4274d537fb..f025d4d132dc 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -5881,20 +5881,25 @@ def _rvs(self, df, size=None, random_state=None):
         return random_state.standard_t(df, size=size)
 
     def _pdf(self, x, df):
-        #                                gamma((df+1)/2)
-        # t.pdf(x, df) = ---------------------------------------------------
-        #                sqrt(pi*df) * gamma(df/2) * (1+x**2/df)**((df+1)/2)
-        r = np.asarray(df*1.0)
-        Px = (np.exp(sc.gammaln((r+1)/2)-sc.gammaln(r/2))
-              / (np.sqrt(r*np.pi)*(1+(x**2)/r)**((r+1)/2)))
-
-        return Px
+        return _lazywhere(
+            df == np.inf, (x, df),
+            f=lambda x, df: norm._pdf(x),
+            f2=lambda x, df: (
+                np.exp(sc.gammaln((df+1)/2)-sc.gammaln(df/2))
+                / (np.sqrt(df*np.pi)*(1+(x**2)/df)**((df+1)/2))
+            )
+        )
 
     def _logpdf(self, x, df):
-        r = df*1.0
-        lPx = (sc.gammaln((r+1)/2) - sc.gammaln(r/2)
-               - (0.5*np.log(r*np.pi) + (r+1)/2*np.log(1+(x**2)/r)))
-        return lPx
+        return _lazywhere(
+            df == np.inf, (x, df),
+            f=lambda x, df: norm._logpdf(x),
+            f2=lambda x, df: (
+                sc.gammaln((df+1)/2) - sc.gammaln(df/2)
+                - (0.5*np.log(df*np.pi)
+                + (df+1)/2*np.log(1+(x**2)/df))
+            )
+        )
 
     def _cdf(self, x, df):
         return sc.stdtr(df, x)
@@ -5909,19 +5914,33 @@ def _isf(self, q, df):
         return -sc.stdtrit(df, q)
 
     def _stats(self, df):
-        mu = np.where(df > 1, 0.0, np.inf)
-        mu2 = _lazywhere(df > 2, (df,),
-                         lambda df: df / (df-2.0),
-                         np.inf)
-        mu2 = np.where(df <= 1, np.nan, mu2)
-        g1 = np.where(df > 3, 0.0, np.nan)
-        g2 = _lazywhere(df > 4, (df,),
-                        lambda df: 6.0 / (df-4.0),
-                        np.inf)
-        g2 = np.where(df <= 2, np.nan, g2)
+        # 0.0, 1.0, 0.0, 0.0
+        cond = df == np.inf
+        mu = np.where(df > 1 | cond, 0.0, np.inf)
+
+        def _mu2_finitedf(df):
+            mu2 = _lazywhere(df > 2, (df,),
+                             lambda df: df / (df-2.0),
+                             np.inf)
+            mu2 = np.where(df <= 1, np.nan, mu2)
+            return mu2
+
+        mu2 = _lazywhere(~cond, (df,), _mu2_finitedf, 1.0)
+        g1 = np.where(df > 3 | cond, 0.0, np.nan)
+
+        def _g2_finitedf(df):
+            g2 = _lazywhere(df > 4, (df,),
+                            lambda df: 6.0 / (df-4.0),
+                            np.inf)
+            g2 = np.where(df <= 2, np.nan, g2)
+            return g2
+
+        g2 = _lazywhere(~cond, (df,), _g2_finitedf, 0.0)
         return mu, mu2, g1, g2
 
     def _entropy(self, df):
+        if df == np.inf:
+            return norm._entropy()
         half = df/2
         half1 = (df + 1)/2
         return (half1*(sc.digamma(half1) - sc.digamma(half))
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 127d65ca22b1..ca033c2efe26 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -2316,6 +2316,48 @@ def test_t_entropy():
     assert_allclose(stats.t.entropy(df), expected, rtol=1e-13)
 
 
+@pytest.mark.parametrize("methname", ["pdf", "logpdf", "cdf",
+                                      "ppf", "sf", "isf"])
+def test_t_inf_df(methname):
+    t_dist = stats.t(df=np.inf, loc=3, scale=1)
+    norm_dist = stats.norm(loc=3, scale=1)
+    t_meth = getattr(t_dist, methname)
+    norm_meth = getattr(norm_dist, methname)
+    assert_equal(t_meth(0.5), norm_meth(0.5))
+    t_dist = stats.t(df=[[np.inf, 1.], [2., np.inf]],
+                     loc=3, scale=1)
+    norm_dist = stats.norm(loc=3, scale=1)
+    t_meth = getattr(t_dist, methname)
+    norm_meth = getattr(norm_dist, methname)
+    res = t_meth([0.5, 0.7])
+    assert_equal(res[0][0], norm_meth(0.5))
+    assert_equal(res[1][1], norm_meth(0.7))
+
+
+def test_t_inf_df_stats_entropy():
+    assert_equal(stats.t.stats(df=np.inf, loc=3, scale=1),
+                 stats.norm.stats(loc=3, scale=1))
+    assert_equal(stats.t.entropy(df=np.inf, loc=3, scale=1),
+                 stats.norm.entropy(loc=3, scale=1))
+    t_dist = stats.t(df=[[np.inf, 1.], [2., np.inf]],
+                     loc=3, scale=1)
+    norm_dist = stats.norm(loc=3, scale=1)
+    res_stats = t_dist.stats('mvsk')
+    res_entropy = t_dist.entropy()
+    assert_equal([res_stats[i][0][0] for i in range(4)],
+                 norm_dist.stats('mvsk'))
+    assert_equal([res_stats[i][1][1] for i in range(4)],
+                 norm_dist.stats('mvsk'))
+    assert_equal([res_stats[i][0][1] for i in range(4)],
+                 stats.t.stats(df=1, loc=3, moments='mvsk'))
+    assert_equal([res_stats[i][1][0] for i in range(4)],
+                 stats.t.stats(df=2, loc=3, moments='mvsk'))
+    assert_equal(res_entropy[0][0], norm_dist.entropy())
+    assert_equal(res_entropy[1][1], norm_dist.entropy())
+    assert_equal(res_entropy[0][1], stats.t.entropy(df=1, loc=3))
+    assert_equal(res_entropy[1][0], stats.t.entropy(df=2, loc=3))
+
+
 class TestRvDiscrete:
     def setup_method(self):
         np.random.seed(1234)


Code review comments:
diff_hunk: @@ -5909,19 +5914,33 @@ def _isf(self, q, df):
         return -sc.stdtrit(df, q)
 
     def _stats(self, df):
-        mu = np.where(df > 1, 0.0, np.inf)
-        mu2 = _lazywhere(df > 2, (df,),
-                         lambda df: df / (df-2.0),
-                         np.inf)
-        mu2 = np.where(df <= 1, np.nan, mu2)
-        g1 = np.where(df > 3, 0.0, np.nan)
-        g2 = _lazywhere(df > 4, (df,),
-                        lambda df: 6.0 / (df-4.0),
-                        np.inf)
-        g2 = np.where(df <= 2, np.nan, g2)
+        # 0.0, 1.0, 0.0, 0.0
+        cond = df == np.inf
+        mu = np.where(df > 1 | cond, 0.0, np.inf)
path: scipy/stats/_continuous_distns.py
text: @user1:
```suggestion
        mu = np.where(df > 1, 0.0, np.inf)
```
`np.inf > 1`, right?

diff_hunk: @@ -5909,19 +5914,33 @@ def _isf(self, q, df):
         return -sc.stdtrit(df, q)
 
     def _stats(self, df):
-        mu = np.where(df > 1, 0.0, np.inf)
-        mu2 = _lazywhere(df > 2, (df,),
-                         lambda df: df / (df-2.0),
-                         np.inf)
-        mu2 = np.where(df <= 1, np.nan, mu2)
-        g1 = np.where(df > 3, 0.0, np.nan)
-        g2 = _lazywhere(df > 4, (df,),
-                        lambda df: 6.0 / (df-4.0),
-                        np.inf)
-        g2 = np.where(df <= 2, np.nan, g2)
+        # 0.0, 1.0, 0.0, 0.0
+        cond = df == np.inf
+        mu = np.where(df > 1 | cond, 0.0, np.inf)
+
+        def _mu2_finitedf(df):
+            mu2 = _lazywhere(df > 2, (df,),
+                             lambda df: df / (df-2.0),
+                             np.inf)
+            mu2 = np.where(df <= 1, np.nan, mu2)
+            return mu2
+
+        mu2 = _lazywhere(~cond, (df,), _mu2_finitedf, 1.0)
path: scipy/stats/_continuous_distns.py
text: @user1:
Please consider tirthasheshpatel/scipy#11.

@author:
Done.