Repository: scipy/scipy
Pull request: ENH: stats.pearsonr: add support for `axis` argument
#### Reference issue
Closes gh-9307

#### What does this implement/fix?
This adds native support for an `axis` argument to `pearsonr`, greatly improving the speed of vectorized calls (e.g. [pairwise comparisons](https://stackoverflow.com/questions/62817623/improve-speed-of-scipy-pearson-correlation-for-many-pairwise-calculations/78025304#78025304)).

#### Additional information
~~Failing existing tests pending resolution of gh-20136.~~ Worked around it.
~~Needs new tests of `axis`.~~ Added.


Problem statement:
feature request: make `scipy.stats.pearsonr` accept 2-D arrays
I'd love to see `scipy.stats.pearsonr` modified to allow 2-D arrays to be passed, with one column per variable and one row per observation.  The output would be a pair of 2-D arrays--one for the Pearson correlation coefficient and one for the p-value.  Currently, obtaining something like this requires a pair of nested loops, with `scipy.stats.pearsonr` being called once for each pair of variables.

Hints from issue discussion:
I know this is probably no more relevant for you, but maybe it could help someone else.
you can convert your data to a pandas dataFrame, and apply the function pandas.DataFrame.corr with parameter method='pearson'. https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.corr.html 
The only (regrettable) difference is that this function does not return p-val. 
I still do such tests.  Thank you!
At some point @swallan or @DominicChm might be vectorizing `mannwhitneyu` to respond to gh-12837; perhaps this would be a good follow-up.
@Phillip-M-Feldman would gh-13312 close this issue, from your perspective? It allows `pearsonr` to process nd-arrays, but it is just for convenience, not speed.
This sounds fine.

On Wed, Dec 30, 2020 at 6:00 PM Matt Haberland <notifications@github.com>
wrote:

> @Phillip-M-Feldman <https://github.com/Phillip-M-Feldman> would gh-13312
> <https://github.com/scipy/scipy/pull/13312> close this issue, from your
> perspective? It allows pearsonr to process nd-arrays, but it is just for
> convenience, not speed.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/scipy/scipy/issues/9307#issuecomment-752815827>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AAIEDRF5EVUKOL5FJBVKYJTSXPLNLANCNFSM4FW7T5RQ>
> .
>

I don't think gh-13312 would address this as requested, actually. I think the request is for it to work more like `spearmanr`, which is an oddball in terms of the way it's vectorized.
On the other hand, if the decorator in gh-13312 were applied to this function, I think the user could achieve the desired behavior by passing in arrays of shapes `(n_var, 1, n_obs)` and `(1, n_var, n_obs)` and taking the statistic along the last axis.

I won't add `pearsonr` in gh-13312, but it could be vectorized like this in a follow-up. 
@tirthasheshpatel @tupui @Kai-Striega thought you might be interested in this one as we try to make the `stats` interface more consistent. We might think that `pearsonr` should be vectorized like `spearmanr`, as they are both correlation functions. But the behavior of `spearmanr` (and `np.corrcoef`) is unusual:

```python3
from scipy.stats import mannwhitneyu, spearmanr
x = np.random.rand(3, 10); y = np.random.rand(3, 10)
mannwhitneyu(x, y, axis=1).statistic.shape  # (3,)
spearmanr(x, y, axis=1)[0].shape  # (6, 6)
```

It does pairwise comparisons between every row of `x` and `y` with every other row of `x` and `y`. 

If we hope to make the `scipy.stats` interface more consistent, I think we'd want it to be vectorized like essentially other function; e.g.
```python3
pearsonr(x, y, axis=1)[0].shape  # we would want (3,), not (6, 6)
```
If the user wants to do pairwise comparisons, then they can do:
```python3
z = np.concatenate([x, y], axis=0)
pearsonr(z[:, np.newaxis, :], z[np.newaxis, :, :], axis=2)[0].shape
```

I'd suggest that we deprecate the 2D behavior then, after the deprecation cycle, replace it.
What do you think? 

@tirthasheshpatel @tupui @Kai-Striega what do you think? Shall we get the deprecation in before 1.9.0 branches so we can vectorize these consistently?
I'm +1 on this but, won't have the bandwidth to implement it

Commit message: MAINT: stats.pearsonr: respect dtype

Patch to review:
diff --git a/scipy/stats/_morestats.py b/scipy/stats/_morestats.py
index 8da934504184..13d6824f7461 100644
--- a/scipy/stats/_morestats.py
+++ b/scipy/stats/_morestats.py
@@ -1250,6 +1250,12 @@ def boxcox_normmax(
     6.000...
     """
     x = np.asarray(x)
+
+    if np.any(~np.isfinite(x) | (x < 0)):
+        message = ("The `x` argument of `boxcox_normmax` must contain "
+                   "only positive, finite, real numbers.")
+        raise ValueError(message)
+
     end_msg = "exceed specified `ymax`."
     if isinstance(ymax, _BigFloat):
         dtype = x.dtype if np.issubdtype(x.dtype, np.floating) else np.float64
@@ -1321,15 +1327,7 @@ def _all(x):
 
     optimfunc = methods[method]
 
-    try:
-        res = optimfunc(x)
-    except ValueError as e:
-        if "infs or NaNs" in str(e):
-            message = ("The `x` argument of `boxcox_normmax` must contain "
-                       "only positive, finite, real numbers.")
-            raise ValueError(message) from e
-        else:
-            raise e
+    res = optimfunc(x)
 
     if res is None:
         message = ("The `optimizer` argument of `boxcox_normmax` must return "
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index a0bcb13c74b2..ead8b5fae5eb 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -45,7 +45,9 @@
                               AxisError)
 
 import scipy.special as special
-from scipy import linalg
+# Import unused here but needs to stay until end of deprecation periode
+# See https://github.com/scipy/scipy/issues/15765#issuecomment-1875564522
+from scipy import linalg  # noqa: F401
 from . import distributions
 from . import _mstats_basic as mstats_basic
 from ._stats_mstats_common import (_find_repeats, linregress, theilslopes,
@@ -4347,13 +4349,11 @@ def _pearsonr_fisher_ci(r, n, confidence_level, alternative):
     Fisher's transformation is used to compute the confidence interval
     (https://en.wikipedia.org/wiki/Fisher_transformation).
     """
-    if r == 1:
-        zr = np.inf
-    elif r == -1:
-        zr = -np.inf
-    else:
+
+    with np.errstate(divide='ignore'):
         zr = np.arctanh(r)
 
+    ones = np.ones_like(r)
     if n > 3:
         se = np.sqrt(1 / (n - 3))
         if alternative == "two-sided":
@@ -4366,28 +4366,28 @@ def _pearsonr_fisher_ci(r, n, confidence_level, alternative):
             h = special.ndtri(confidence_level)
             zhi = zr + h*se
             rhi = np.tanh(zhi)
-            rlo = -1.0
+            rlo = -ones
         else:
             # alternative == "greater":
             h = special.ndtri(confidence_level)
             zlo = zr - h*se
             rlo = np.tanh(zlo)
-            rhi = 1.0
+            rhi = ones
     else:
-        rlo, rhi = -1.0, 1.0
+        rlo, rhi = -ones, ones
 
     return ConfidenceInterval(low=rlo, high=rhi)
 
 
-def _pearsonr_bootstrap_ci(confidence_level, method, x, y, alternative):
+def _pearsonr_bootstrap_ci(confidence_level, method, x, y, alternative, axis):
     """
     Compute the confidence interval for Pearson's R using the bootstrap.
     """
-    def statistic(x, y):
-        statistic, _ = pearsonr(x, y)
+    def statistic(x, y, axis):
+        statistic, _ = pearsonr(x, y, axis=axis)
         return statistic
 
-    res = bootstrap((x, y), statistic, confidence_level=confidence_level,
+    res = bootstrap((x, y), statistic, confidence_level=confidence_level, axis=axis,
                     paired=True, alternative=alternative, **method._asdict())
     # for one-sided confidence intervals, bootstrap gives +/- inf on one side
     res.confidence_interval = np.clip(res.confidence_interval, -1, 1)
@@ -4419,12 +4419,13 @@ class PearsonRResult(PearsonRResultBase):
         coefficient `statistic` for the given confidence level.
 
     """
-    def __init__(self, statistic, pvalue, alternative, n, x, y):
+    def __init__(self, statistic, pvalue, alternative, n, x, y, axis):
         super().__init__(statistic, pvalue)
         self._alternative = alternative
         self._n = n
         self._x = x
         self._y = y
+        self._axis = axis
 
         # add alias for consistency with other correlation functions
         self.correlation = statistic
@@ -4475,8 +4476,8 @@ def confidence_interval(self, confidence_level=0.95, method=None):
                https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
         """
         if isinstance(method, BootstrapMethod):
-            ci = _pearsonr_bootstrap_ci(confidence_level, method,
-                                        self._x, self._y, self._alternative)
+            ci = _pearsonr_bootstrap_ci(confidence_level, method, self._x, self._y,
+                                        self._alternative, self._axis)
         elif method is None:
             ci = _pearsonr_fisher_ci(self.statistic, self._n, confidence_level,
                                      self._alternative)
@@ -4486,7 +4487,7 @@ def confidence_interval(self, confidence_level=0.95, method=None):
             raise ValueError(message)
         return ci
 
-def pearsonr(x, y, *, alternative='two-sided', method=None):
+def pearsonr(x, y, *, alternative='two-sided', method=None, axis=0):
     r"""
     Pearson correlation coefficient and p-value for testing non-correlation.
 
@@ -4508,10 +4509,15 @@ def pearsonr(x, y, *, alternative='two-sided', method=None):
 
     Parameters
     ----------
-    x : (N,) array_like
+    x : array_like
         Input array.
-    y : (N,) array_like
+    y : array_like
         Input array.
+    axis : int or None, default
+        Axis along which to perform the calculation. Default is 0.
+        If None, ravel both arrays before performing the calculation.
+
+        .. versionadded:: 1.13.0
     alternative : {'two-sided', 'greater', 'less'}, optional
         Defines the alternative hypothesis. Default is 'two-sided'.
         The following options are available:
@@ -4669,6 +4675,24 @@ def pearsonr(x, y, *, alternative='two-sided', method=None):
     >>> res.confidence_interval(confidence_level=0.9, method=method)
     ConfidenceInterval(low=-0.9983163756488651, high=-0.22771001702132443)  # may vary
 
+    If N-dimensional arrays are provided, multiple tests are performed in a
+    single call according to the same conventions as most `scipy.stats` functions:
+
+    >>> rng = np.random.default_rng(2348246935601934321)
+    >>> x = rng.standard_normal((8, 15))
+    >>> y = rng.standard_normal((8, 15))
+    >>> stats.pearsonr(x, y, axis=0).statistic.shape  # between corresponding columns
+    (15,)
+    >>> stats.pearsonr(x, y, axis=1).statistic.shape  # between corresponding rows
+    (8,)
+
+    To perform all pairwise comparisons between slices of the arrays,
+    use standard NumPy broadcasting techniques. For instance, to compute the
+    correlation between all pairs of rows:
+
+    >>> stats.pearsonr(x[:, np.newaxis, :], y, axis=-1).statistic.shape
+    (8, 8)
+
     There is a linear dependence between x and y if y = a + b*x + e, where
     a,b are constants and e is a random error term, assumed to be independent
     of x. For simplicity, assume that x is standard normal, a=0, b=1 and let
@@ -4716,85 +4740,115 @@ def pearsonr(x, y, *, alternative='two-sided', method=None):
     than zero which happens in about half of the cases if we sample x and y.
 
     """
-    n = len(x)
-    if n != len(y):
-        raise ValueError('x and y must have the same length.')
+    x = np.asarray(x)
+    y = np.asarray(y)
+
+    if axis is None:
+        x = x.ravel()
+        y = y.ravel()
+        axis = -1
+
+    axis_int = int(axis)
+    if axis_int != axis:
+        raise ValueError('`axis` must be an integer.')
+    axis = axis_int
+
+    n = x.shape[axis]
+    if n != y.shape[axis]:
+        raise ValueError('`x` and `y` must have the same length along `axis`.')
 
     if n < 2:
-        raise ValueError('x and y must have length at least 2.')
+        raise ValueError('`x` and `y` must have length at least 2.')
 
-    x = np.asarray(x)
-    y = np.asarray(y)
+    try:
+        x, y = np.broadcast_arrays(x, y)
+    except ValueError as e:
+        message = '`x` and `y` must be broadcastable.'
+        raise ValueError(message) from e
 
-    if (np.issubdtype(x.dtype, np.complexfloating)
-            or np.issubdtype(y.dtype, np.complexfloating)):
+    x = np.moveaxis(x, axis, -1)
+    y = np.moveaxis(y, axis, -1)
+    axis = -1
+
+    dtype = np.result_type(x.dtype, y.dtype)
+    if np.issubdtype(dtype, np.integer):
+        dtype = np.float64
+
+    if np.issubdtype(dtype, np.complexfloating):
         raise ValueError('This function does not support complex data')
 
+    x = x.astype(np.float64, copy=False)
+    y = y.astype(np.float64, copy=False)
+    threshold = np.finfo(dtype).eps ** 0.75
+
     # If an input is constant, the correlation coefficient is not defined.
-    if (x == x[0]).all() or (y == y[0]).all():
+    const_x = (x == x[..., [0]]).all(axis=-1)
+    const_y = (y == y[..., [0]]).all(axis=-1)
+    const_xy = const_x | const_y
+    if np.any(const_xy):
         msg = ("An input array is constant; the correlation coefficient "
                "is not defined.")
         warnings.warn(stats.ConstantInputWarning(msg), stacklevel=2)
-        result = PearsonRResult(statistic=np.nan, pvalue=np.nan, n=n,
-                                alternative=alternative, x=x, y=y)
-        return result
 
     if isinstance(method, PermutationMethod):
-        def statistic(y):
-            statistic, _ = pearsonr(x, y, alternative=alternative)
+        def statistic(y, axis):
+            statistic, _ = pearsonr(x, y, axis=axis, alternative=alternative)
             return statistic
 
         res = permutation_test((y,), statistic, permutation_type='pairings',
-                               alternative=alternative, **method._asdict())
+                               axis=axis, alternative=alternative, **method._asdict())
 
         return PearsonRResult(statistic=res.statistic, pvalue=res.pvalue, n=n,
-                              alternative=alternative, x=x, y=y)
+                              alternative=alternative, x=x, y=y, axis=axis)
     elif isinstance(method, MonteCarloMethod):
-        def statistic(x, y):
-            statistic, _ = pearsonr(x, y, alternative=alternative)
+        def statistic(x, y, axis):
+            statistic, _ = pearsonr(x, y, axis=axis, alternative=alternative)
             return statistic
 
         if method.rvs is None:
             rng = np.random.default_rng()
             method.rvs = rng.normal, rng.normal
 
-        res = monte_carlo_test((x, y,), statistic=statistic,
+        res = monte_carlo_test((x, y,), statistic=statistic, axis=axis,
                                alternative=alternative, **method._asdict())
 
         return PearsonRResult(statistic=res.statistic, pvalue=res.pvalue, n=n,
-                              alternative=alternative, x=x, y=y)
+                              alternative=alternative, x=x, y=y, axis=axis)
     elif method is not None:
         message = ('`method` must be an instance of `PermutationMethod`,'
                    '`MonteCarloMethod`, or None.')
         raise ValueError(message)
 
-    # dtype is the data type for the calculations.  This expression ensures
-    # that the data type is at least 64 bit floating point.  It might have
-    # more precision if the input is, for example, np.longdouble.
-    dtype = type(1.0 + x[0] + y[0])
-
     if n == 2:
-        r = dtype(np.sign(x[1] - x[0])*np.sign(y[1] - y[0]))
-        result = PearsonRResult(statistic=r, pvalue=1.0, n=n,
-                                alternative=alternative, x=x, y=y)
+        r = np.asarray(np.sign(x[..., 1] - x[..., 0])*np.sign(y[..., 1] - y[..., 0]),
+                       dtype=dtype)
+        result = PearsonRResult(statistic=r, pvalue=np.ones_like(r)[()], n=n,
+                                alternative=alternative, x=x, y=y, axis=axis)
         return result
 
-    xmean = x.mean(dtype=dtype)
-    ymean = y.mean(dtype=dtype)
+    xmean = x.mean(axis=axis, keepdims=True, dtype=dtype)
+    ymean = y.mean(axis=axis, keepdims=True, dtype=dtype)
 
     # By using `astype(dtype)`, we ensure that the intermediate calculations
     # use at least 64 bit floating point.
     xm = x.astype(dtype) - xmean
     ym = y.astype(dtype) - ymean
 
-    # Unlike np.linalg.norm or the expression sqrt((xm*xm).sum()),
-    # scipy.linalg.norm(xm) does not overflow if xm is, for example,
+    # scipy.linalg.norm(xm) avoids premature overflow when xm is e.g.
     # [-5e210, 5e210, 3e200, -3e200]
-    normxm = linalg.norm(xm)
-    normym = linalg.norm(ym)
-
-    threshold = 1e-13
-    if normxm < threshold*abs(xmean) or normym < threshold*abs(ymean):
+    # but not when `axis` is provided, so scale manually. scipy.linalg.norm
+    # also raises an error with NaN input rather than returning NaN, so
+    # use np.linalg.norm.
+    xmax = np.max(np.abs(xm), axis=axis, keepdims=True)
+    ymax = np.max(np.abs(ym), axis=axis, keepdims=True)
+    with np.errstate(invalid='ignore'):
+        normxm = xmax * np.linalg.norm(xm/xmax, axis=axis, keepdims=True)
+        normym = ymax * np.linalg.norm(ym/ymax, axis=axis, keepdims=True)
+
+    nconst_x = np.any(normxm < threshold*abs(xmean), axis=axis)
+    nconst_y = np.any(normym < threshold*abs(ymean), axis=axis)
+    nconst_xy = nconst_x | nconst_y
+    if np.any(nconst_xy & (~const_xy)):
         # If all the values in x (likewise y) are very close to the mean,
         # the loss of precision that occurs in the subtraction xm = x - xmean
         # might result in large errors in r.
@@ -4802,20 +4856,22 @@ def statistic(x, y):
                "correlation coefficient may be inaccurate.")
         warnings.warn(stats.NearConstantInputWarning(msg), stacklevel=2)
 
-    r = np.dot(xm/normxm, ym/normym)
+    with np.errstate(invalid='ignore', divide='ignore'):
+        r = np.sum(xm/normxm * ym/normym, axis=axis)
 
     # Presumably, if abs(r) > 1, then it is only some small artifact of
     # floating point arithmetic.
-    r = max(min(r, 1.0), -1.0)
+    r = np.asarray(np.clip(r, -1, 1))
+    r[const_xy] = np.nan
 
     # As explained in the docstring, the distribution of `r` under the null
     # hypothesis is the beta distribution on (-1, 1) with a = b = n/2 - 1.
-    ab = n/2 - 1
-    dist = stats.beta(ab, ab, loc=-1, scale=2)
+    ab, loc, scale = np.asarray([n/2 - 1, -1, 2], dtype=dtype)
+    dist = stats.beta(ab, ab, loc=loc, scale=scale)
     pvalue = _get_pvalue(r, dist, alternative)
 
-    return PearsonRResult(statistic=r, pvalue=pvalue, n=n,
-                          alternative=alternative, x=x, y=y)
+    return PearsonRResult(statistic=r, pvalue=pvalue.astype(dtype), n=n,
+                          alternative=alternative, x=x, y=y, axis=axis)
 
 
 def fisher_exact(table, alternative='two-sided'):
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index a1600a95e560..ab7e2ed71e0d 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -401,7 +401,7 @@ def test_very_small_input_values(self):
         # See https://github.com/scipy/scipy/issues/9353
         x = [0.004434375, 0.004756007, 0.003911996, 0.0038005, 0.003409971]
         y = [2.48e-188, 7.41e-181, 4.09e-208, 2.08e-223, 2.66e-245]
-        r, p = stats.pearsonr(x,y)
+        r, p = stats.pearsonr(x, y)
 
         # The expected values were computed using mpmath with 80 digits
         # of precision.
@@ -516,30 +516,30 @@ def test_complex_data(self):
     @pytest.mark.parametrize('method', ('permutation', 'monte_carlo'))
     def test_resampling_pvalue(self, method, alternative):
         rng = np.random.default_rng(24623935790378923)
-        size = 100 if method == 'permutation' else 1000
+        size = (2, 100) if method == 'permutation' else (2, 1000)
         x = rng.normal(size=size)
         y = rng.normal(size=size)
         methods = {'permutation': stats.PermutationMethod(random_state=rng),
                    'monte_carlo': stats.MonteCarloMethod(rvs=(rng.normal,)*2)}
         method = methods[method]
-        res = stats.pearsonr(x, y, alternative=alternative, method=method)
-        ref = stats.pearsonr(x, y, alternative=alternative)
+        res = stats.pearsonr(x, y, alternative=alternative, method=method, axis=-1)
+        ref = stats.pearsonr(x, y, alternative=alternative, axis=-1)
         assert_allclose(res.statistic, ref.statistic, rtol=1e-15)
         assert_allclose(res.pvalue, ref.pvalue, rtol=1e-2, atol=1e-3)
 
     @pytest.mark.xslow
     @pytest.mark.parametrize('alternative', ('less', 'greater', 'two-sided'))
     def test_bootstrap_ci(self, alternative):
-        rng = np.random.default_rng(24623935790378923)
-        x = rng.normal(size=100)
-        y = rng.normal(size=100)
-        res = stats.pearsonr(x, y, alternative=alternative)
+        rng = np.random.default_rng(2462935790378923)
+        x = rng.normal(size=(2, 100))
+        y = rng.normal(size=(2, 100))
+        res = stats.pearsonr(x, y, alternative=alternative, axis=-1)
 
         method = stats.BootstrapMethod(random_state=rng)
         res_ci = res.confidence_interval(method=method)
         ref_ci = res.confidence_interval()
 
-        assert_allclose(res_ci, ref_ci, atol=1e-2)
+        assert_allclose(res_ci, ref_ci, atol=1.5e-2)
 
     def test_invalid_method(self):
         message = "`method` must be an instance of..."
@@ -550,6 +550,89 @@ def test_invalid_method(self):
         with pytest.raises(ValueError, match=message):
             res.confidence_interval(method="exact")
 
+    @pytest.mark.parametrize('axis', [0, 1])
+    def test_axis01(self, axis):
+        rng = np.random.default_rng(38572345825)
+        shape = (9, 10)
+        x, y = rng.normal(size=(2,) + shape)
+        res = stats.pearsonr(x, y, axis=axis)
+        ci = res.confidence_interval()
+        if axis == 0:
+            x, y = x.T, y.T
+        for i in range(x.shape[0]):
+            res_i = stats.pearsonr(x[i], y[i])
+            ci_i = res_i.confidence_interval()
+            assert_allclose(res.statistic[i], res_i.statistic)
+            assert_allclose(res.pvalue[i], res_i.pvalue)
+            assert_allclose(ci.low[i], ci_i.low)
+            assert_allclose(ci.high[i], ci_i.high)
+
+    def test_axis_None(self):
+        rng = np.random.default_rng(38572345825)
+        shape = (9, 10)
+        x, y = rng.normal(size=(2,) + shape)
+        res = stats.pearsonr(x, y, axis=None)
+        ci = res.confidence_interval()
+        ref = stats.pearsonr(x.ravel(), y.ravel())
+        ci_ref = ref.confidence_interval()
+        assert_allclose(res.statistic, ref.statistic)
+        assert_allclose(res.pvalue, ref.pvalue)
+        assert_allclose(ci, ci_ref)
+
+    def test_nd_input_validation(self):
+        x, y = np.ones((2, 2, 5))
+        message = '`axis` must be an integer.'
+        with pytest.raises(ValueError, match=message):
+            stats.pearsonr(x, y, axis=1.5)
+
+        message = '`x` and `y` must have the same length along `axis`'
+        with pytest.raises(ValueError, match=message):
+            stats.pearsonr(x, np.ones((2, 1)), axis=1)
+
+        message = '`x` and `y` must have length at least 2.'
+        with pytest.raises(ValueError, match=message):
+            stats.pearsonr(*np.ones((2, 2, 1)), axis=1)
+
+        message = '`x` and `y` must be broadcastable.'
+        with pytest.raises(ValueError, match=message):
+            stats.pearsonr(x, np.ones((3, 5)), axis=1)
+
+    def test_nd_special_cases(self):
+        rng = np.random.default_rng(34989235492245)
+        x0, y0 = rng.random((2, 3, 5))
+
+        message = 'An input array is constant'
+        with pytest.warns(stats.ConstantInputWarning, match=message):
+            x = x0.copy()
+            x[0] = 1
+            res = stats.pearsonr(x, y0, axis=1)
+            ci = res.confidence_interval()
+            np.testing.assert_equal(res.statistic[0], np.nan)
+            np.testing.assert_equal(res.pvalue[0], np.nan)
+            np.testing.assert_equal(ci.low[0], np.nan)
+            np.testing.assert_equal(ci.high[0], np.nan)
+            assert not np.any(np.isnan(res.statistic[1:]))
+            assert not np.any(np.isnan(res.pvalue[1:]))
+            assert not np.any(np.isnan(ci.low[1:]))
+            assert not np.any(np.isnan(ci.high[1:]))
+
+        message = 'An input array is nearly constant'
+        with pytest.warns(stats.NearConstantInputWarning, match=message):
+            x[0, 0] = 1 + 1e-15
+            stats.pearsonr(x, y0, axis=1)
+
+        # length 2 along axis
+        x = [[1, 2], [1, 2], [2, 1], [2, 1]]
+        y = [[1, 2], [2, 1], [1, 2], [2, 1]]
+        res = stats.pearsonr(x, y, axis=-1)
+        ci = res.confidence_interval()
+        assert_allclose(res.statistic, [1, -1, -1, 1])
+        assert_allclose(res.pvalue, 1)
+        assert_allclose(ci.low, -1)
+        assert_allclose(ci.high, 1)
+        assert (res.pvalue.shape == res.statistic.shape
+                == ci.low.shape == ci.high.shape)
+
 
 class TestFisherExact:
     """Some tests to show that fisher_exact() works correctly.


Code review comments:
diff_hunk: @@ -4716,106 +4740,138 @@ def pearsonr(x, y, *, alternative='two-sided', method=None):
     than zero which happens in about half of the cases if we sample x and y.
 
     """
-    n = len(x)
-    if n != len(y):
-        raise ValueError('x and y must have the same length.')
+    x = np.asarray(x)
+    y = np.asarray(y)
+
+    if axis is None:
+        x = x.ravel()
+        y = y.ravel()
+        axis = -1
+
+    axis_int = int(axis)
+    if axis_int != axis:
+        raise ValueError('`axis` must be an integer.')
+    axis = axis_int
+
+    n = x.shape[axis]
+    if n != y.shape[axis]:
+        raise ValueError('`x` and `y` must have the same length along `axis`.')
 
     if n < 2:
-        raise ValueError('x and y must have length at least 2.')
+        raise ValueError('`x` and `y` must have length at least 2.')
 
-    x = np.asarray(x)
-    y = np.asarray(y)
+    try:
+        x, y = np.broadcast_arrays(x, y)
+    except ValueError as e:
+        message = '`x` and `y` must be broadcastable.'
+        raise ValueError(message) from e
 
-    if (np.issubdtype(x.dtype, np.complexfloating)
-            or np.issubdtype(y.dtype, np.complexfloating)):
+    x = np.moveaxis(x, axis, -1)
+    y = np.moveaxis(y, axis, -1)
+    axis = -1
+
+    dtype = np.result_type(x.dtype, y.dtype)
+    if np.issubdtype(dtype, np.integer):
+        dtype = np.float64
+
+    if np.issubdtype(dtype, np.complexfloating):
         raise ValueError('This function does not support complex data')
 
+    x = x.astype(np.float64, copy=False)
+    y = y.astype(np.float64, copy=False)
path: scipy/stats/_stats_py.py
text: @author:
```suggestion
    x = x.astype(dtype, copy=False)
    y = y.astype(dtype, copy=False)
```