Repository: scipy/scipy
Pull request: MAINT: stats.sampling.NumericalInverseHermite: private distribution methods need at least 1d arrays
#### Reference issue
Closes gh-17155

#### What does this implement/fix?
This fixes gh-17155.

#### Additional information
I'm not sure if there are other implications of this change. Based on the code like:
```python3
   if res < 0:
       return 0
```
It looks like `res` is _always_ expected to have one element. So since I'm ensuring that it will be 1d with one element, shall I extract that one element from the result?


Problem statement:
BUG: stats: Bug in XSLOW tests in TestNumericalInverseHermite
### Describe your issue.

When the environment variable `SCIPY_XSLOW` is set, I get two failures of a test in the `TestNumericalInverseHermite` test class.  I am also use a development version of numpy (version `'1.24.0.dev0+837.g9c5a2b41c'`).

Here's the transcript of running just the `TestNumericalInverseHermite` tests:

```
$ SCIPY_XSLOW=1 python dev.py test -m full -j 8 -t scipy.stats.tests.test_sampling::TestNumericalInverseHermite
💻  ninja -C /home/warren/repos/git/forks/scipy/build
ninja: Entering directory `/home/warren/repos/git/forks/scipy/build'
[2/2] Generating scipy/generate-config with a custom command
Build OK
💻  meson install -C build --only-changed
Installing, see meson-install.log...
Installation OK
SciPy from development installed path at: /home/warren/repos/git/forks/scipy/build-install/lib/python3.10/site-packages
Running tests for scipy version:1.10.0.dev0+2006.527f2dd, installed at:/home/warren/repos/git/forks/scipy/build-install/lib/python3.10/site-packages/scipy
============================================================== test session starts ===============================================================
platform linux -- Python 3.10.1, pytest-6.2.5, py-1.11.0, pluggy-1.0.0
rootdir: /home/warren/repos/git/forks/scipy, configfile: pytest.ini
plugins: forked-1.4.0, timeout-2.1.0, hypothesis-6.35.0, xdist-2.5.0
gw0 [145] / gw1 [145] / gw2 [145] / gw3 [145] / gw4 [145] / gw5 [145] / gw6 [145] / gw7 [145]
...............x............x.......x.........xxx....ss.x...x..............x.x...........xx.xs.s..........x.........................F....x [ 95%]
..F....                                                                                                                                    [100%]
==================================================================== FAILURES ====================================================================
__________________________________ TestNumericalInverseHermite.test_basic_all_scipy_dists[truncnorm-shapes101] ___________________________________
[gw7] linux -- Python 3.10.1 /home/warren/py3.10.1/bin/python
scipy/stats/tests/test_sampling.py:1019: in test_basic_all_scipy_dists
    fni = NumericalInverseHermite(dist)
        dist       = <scipy.stats._distn_infrastructure.rv_continuous_frozen object at 0x7f9e0d8fbc70>
        distname   = 'truncnorm'
        fail_dists = {'beta', 'gausshyper', 'genhyperbolic', 'geninvgauss', 'invgauss', 'kappa4', ...}
        self       = <scipy.stats.tests.test_sampling.TestNumericalInverseHermite object at 0x7f9e0d8fb1c0>
        shapes     = (0.1, 2.0)
        slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
        sup        = <numpy.testing._private.utils.suppress_warnings object at 0x7f9e0d8fbaf0>
scipy/stats/unuran_wrapper.pyx:2020: in scipy.stats.unuran_wrapper.NumericalInverseHermite.__cinit__
    ???
scipy/stats/unuran_wrapper.pyx:1996: in scipy.stats.unuran_wrapper.NumericalInverseHermite.__cinit__._callback_wrapper
    ???
scipy/stats/unuran_wrapper.pyx:271: in scipy.stats.unuran_wrapper._unpack_dist.wrap_dist.cdf
    ???
scipy/stats/_continuous_distns.py:8670: in _cdf
    return np.exp(self._logcdf(x, a, b))
        a          = 0.1
        b          = 2.0
        self       = <scipy.stats._continuous_distns.truncnorm_gen object at 0x7f9e0d8fb190>
        x          = array(2.)
scipy/stats/_continuous_distns.py:8677: in _logcdf
    logcdf[i] = np.log1p(-np.exp(self._logsf(x[i], a[i], b[i])))
E   IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed
        a          = array(0.1)
        b          = array(2.)
        i          = array([ True])
        logcdf     = array([0.])
        self       = <scipy.stats._continuous_distns.truncnorm_gen object at 0x7f9e0d8fb190>
        x          = array(2.)
__________________________________ TestNumericalInverseHermite.test_basic_all_scipy_dists[truncnorm-shapes100] ___________________________________
[gw1] linux -- Python 3.10.1 /home/warren/py3.10.1/bin/python
scipy/stats/tests/test_sampling.py:1019: in test_basic_all_scipy_dists
    fni = NumericalInverseHermite(dist)
        dist       = <scipy.stats._distn_infrastructure.rv_continuous_frozen object at 0x7fb3dff6fd00>
        distname   = 'truncnorm'
        fail_dists = {'beta', 'gausshyper', 'genhyperbolic', 'geninvgauss', 'invgauss', 'kappa4', ...}
        self       = <scipy.stats.tests.test_sampling.TestNumericalInverseHermite object at 0x7fb3dff6f6a0>
        shapes     = (-1.0978730080013919, 2.730675410903198)
        slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
        sup        = <numpy.testing._private.utils.suppress_warnings object at 0x7fb3dff6dd80>
scipy/stats/unuran_wrapper.pyx:2020: in scipy.stats.unuran_wrapper.NumericalInverseHermite.__cinit__
    ???
scipy/stats/unuran_wrapper.pyx:1996: in scipy.stats.unuran_wrapper.NumericalInverseHermite.__cinit__._callback_wrapper
    ???
scipy/stats/unuran_wrapper.pyx:271: in scipy.stats.unuran_wrapper._unpack_dist.wrap_dist.cdf
    ???
scipy/stats/_continuous_distns.py:8670: in _cdf
    return np.exp(self._logcdf(x, a, b))
        a          = -1.0978730080013919
        b          = 2.730675410903198
        self       = <scipy.stats._continuous_distns.truncnorm_gen object at 0x7fb3dff6d7b0>
        x          = array(2.73067541)
scipy/stats/_continuous_distns.py:8677: in _logcdf
    logcdf[i] = np.log1p(-np.exp(self._logsf(x[i], a[i], b[i])))
E   IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed
        a          = array(-1.09787301)
        b          = array(2.73067541)
        i          = array([ True])
        logcdf     = array([0.])
        self       = <scipy.stats._continuous_distns.truncnorm_gen object at 0x7fb3dff6d7b0>
        x          = array(2.73067541)
============================================================ short test summary info =============================================================
FAILED scipy/stats/tests/test_sampling.py::TestNumericalInverseHermite::test_basic_all_scipy_dists[truncnorm-shapes101] - IndexError: too many ...
FAILED scipy/stats/tests/test_sampling.py::TestNumericalInverseHermite::test_basic_all_scipy_dists[truncnorm-shapes100] - IndexError: too many ...
============================================== 2 failed, 124 passed, 4 skipped, 15 xfailed in 3.32s ==============================================
warren@pop-os:~/repos/git/forks/scipy$ 

```

### Reproducing Code Example

```python
n/a
```


### Error message

```shell
n/a
```


### SciPy/NumPy/Python version information

1.10.0.dev0+2006.527f2dd 1.24.0.dev0+837.g9c5a2b41c sys.version_info(major=3, minor=10, micro=1, releaselevel='final', serial=0)

Hints from issue discussion:
The existing distribution infrastructure usually ensures that the inputs to private (distribution-specific implementation) methods are at least 1d, but it is not perfectly consistent. This could be one of those cases, or perhaps NumericalInverseHermite is calling `_logcdf` directly and passing a 0d array? Fix looks simple; I'll give it a shot today.

Commit message: Merge branch 'main' into gh17155

Patch to review:
diff --git a/scipy/stats/_unuran/unuran_wrapper.pyx.templ b/scipy/stats/_unuran/unuran_wrapper.pyx.templ
index 1c44d681b273..90c1a79b3e32 100644
--- a/scipy/stats/_unuran/unuran_wrapper.pyx.templ
+++ b/scipy/stats/_unuran/unuran_wrapper.pyx.templ
@@ -264,10 +264,10 @@ cdef dict _unpack_dist(object dist, str dist_type, list meths,
                     self.support = dist.support
                 def pdf(self, x):
                     # some distributions require array inputs.
-                    x = np.asarray((x-self.loc)/self.scale)
+                    x = np.atleast_1d((x-self.loc)/self.scale)
                     return max(0, self.dist.dist._pdf(x, *self.args)/self.scale)
                 def cdf(self, x):
-                    x = np.asarray((x-self.loc)/self.scale)
+                    x = np.atleast_1d((x-self.loc)/self.scale)
                     res = self.dist.dist._cdf(x, *self.args)
                     if res < 0:
                         return 0
@@ -284,10 +284,10 @@ cdef dict _unpack_dist(object dist, str dist_type, list meths,
                     self.support = dist.support
                 def pmf(self, x):
                     # some distributions require array inputs.
-                    x = np.asarray(x-self.loc)
+                    x = np.atleast_1d(x-self.loc)
                     return max(0, self.dist.dist._pmf(x, *self.args))
                 def cdf(self, x):
-                    x = np.asarray(x-self.loc)
+                    x = np.atleast_1d(x-self.loc)
                     res = self.dist.dist._cdf(x, *self.args)
                     if res < 0:
                         return 0
@@ -797,7 +797,7 @@ cdef class TransformedDensityRejection(Method):
     ...         return 1-x*x
     ...     def dpdf(self, x):
     ...         return -2*x
-    ... 
+    ...
     >>> dist = MyDist()
     >>> rng = TransformedDensityRejection(dist, domain=(-1, 1),
     ...                                   random_state=urng)
@@ -813,7 +813,7 @@ cdef class TransformedDensityRejection(Method):
     ...         return -2*x
     ...     def support(self):
     ...         return (-1, 1)
-    ... 
+    ...
     >>> dist = MyDist()
     >>> rng = TransformedDensityRejection(dist, random_state=urng)
 
@@ -980,16 +980,16 @@ cdef class TransformedDensityRejection(Method):
         >>> from scipy.stats.sampling import TransformedDensityRejection
         >>> from scipy.stats import norm
         >>> from math import exp
-        >>> 
+        >>>
         >>> class MyDist:
         ...     def pdf(self, x):
         ...         return exp(-0.5 * x**2)
         ...     def dpdf(self, x):
         ...         return -x * exp(-0.5 * x**2)
-        ... 
+        ...
         >>> dist = MyDist()
         >>> rng = TransformedDensityRejection(dist)
-        >>> 
+        >>>
         >>> rng.ppf_hat(0.5)
         -0.00018050266342393984
         >>> norm.ppf(0.5)
@@ -1250,7 +1250,7 @@ cdef class NumericalInversePolynomial(Method):
     * The algorithm has problems when the distribution has heavy tails (as then the
       inverse CDF becomes very steep at 0 or 1) and the requested u-resolution is
       very small. E.g., the Cauchy distribution is likely to show this problem when
-      the requested u-resolution is less then 1.e-12. 
+      the requested u-resolution is less then 1.e-12.
 
 
     Parameters
@@ -1321,7 +1321,7 @@ cdef class NumericalInversePolynomial(Method):
     >>> class StandardNormal:
     ...    def pdf(self, x):
     ...        return np.exp(-0.5 * x*x)
-    ... 
+    ...
     >>> dist = StandardNormal()
     >>> urng = np.random.default_rng()
     >>> rng = NumericalInversePolynomial(dist, random_state=urng)
@@ -1357,7 +1357,7 @@ cdef class NumericalInversePolynomial(Method):
     ...        return np.exp(-0.5 * x*x)
     ...    def cdf(self, x):
     ...        return ndtr(x)
-    ... 
+    ...
     >>> dist = StandardNormal()
     >>> urng = np.random.default_rng()
     >>> rng = NumericalInversePolynomial(dist, random_state=urng)
@@ -1721,7 +1721,7 @@ cdef class NumericalInversePolynomial(Method):
             if d == 1:
                 return qrvs.reshape(tuple_size)
             else:
-                return qrvs.reshape(tuple_size + (d,))    
+                return qrvs.reshape(tuple_size + (d,))
 
 
 cdef class NumericalInverseHermite(Method):
@@ -1852,7 +1852,7 @@ cdef class NumericalInverseHermite(Method):
     ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
     ...     def cdf(self, x):
     ...        return ndtr(x)
-    ... 
+    ...
     >>> dist = StandardNormal()
     >>> urng = np.random.default_rng()
     >>> rng = NumericalInverseHermite(dist, random_state=urng)
@@ -1906,7 +1906,7 @@ cdef class NumericalInverseHermite(Method):
     ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
     ...     def cdf(self, x):
     ...        return ndtr(x)
-    ... 
+    ...
     >>> dist = StandardNormal()
     >>> urng = np.random.default_rng()
     >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
@@ -2347,7 +2347,7 @@ cdef class DiscreteAliasUrn(Method):
     ...         return self.p**x * (1-self.p)**(self.n-x)
     ...     def support(self):
     ...         return (0, self.n)
-    ... 
+    ...
     >>> n, p = 10, 0.2
     >>> dist = Binomial(n, p)
     >>> rng = DiscreteAliasUrn(dist, random_state=urng)
diff --git a/scipy/stats/tests/test_sampling.py b/scipy/stats/tests/test_sampling.py
index 1d24dbb1a9c2..c2bd134a4523 100644
--- a/scipy/stats/tests/test_sampling.py
+++ b/scipy/stats/tests/test_sampling.py
@@ -1005,10 +1005,8 @@ def test_inf_nan_domains(self, domain, err, msg):
         with pytest.raises(err, match=msg):
             NumericalInverseHermite(StandardNormal(), domain=domain)
 
-    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
-    @pytest.mark.xslow
-    @pytest.mark.parametrize(("distname", "shapes"), distcont)
-    def test_basic_all_scipy_dists(self, distname, shapes):
+
+    def basic_test_all_scipy_dists(self, distname, shapes):
         slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
         fail_dists = {'beta', 'gausshyper', 'geninvgauss', 'ncf', 'nct',
                       'norminvgauss', 'genhyperbolic', 'studentized_range',
@@ -1033,6 +1031,18 @@ def test_basic_all_scipy_dists(self, distname, shapes):
         assert p_tol < 1e-8
         assert u_tol < 1e-12
 
+    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
+    @pytest.mark.xslow
+    @pytest.mark.parametrize(("distname", "shapes"), distcont)
+    def test_basic_all_scipy_dists(self, distname, shapes):
+        # if distname == "truncnorm":
+        #     pytest.skip("Tested separately")
+        self.basic_test_all_scipy_dists(distname, shapes)
+
+    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
+    def test_basic_truncnorm_gh17155(self):
+        self.basic_test_all_scipy_dists("truncnorm", (0.1, 2))
+
     def test_input_validation(self):
         match = r"`order` must be either 1, 3, or 5."
         with pytest.raises(ValueError, match=match):


Code review comments:
diff_hunk: @@ -1005,10 +1005,8 @@ def test_inf_nan_domains(self, domain, err, msg):
         with pytest.raises(err, match=msg):
             NumericalInverseHermite(StandardNormal(), domain=domain)
 
-    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
-    @pytest.mark.xslow
-    @pytest.mark.parametrize(("distname", "shapes"), distcont)
-    def test_basic_all_scipy_dists(self, distname, shapes):
+
path: scipy/stats/tests/test_sampling.py
text: @author:
```suggestion
```

@author:
If other tests pass, this can probably get committed with [skip ci] before merge to fix the lint issue.

@user1:
Just a small change, so I applied it. Hope you don't mind!

@author:
Yup, that was the intent. Thanks!