Repository: scipy/scipy
Pull request: MAINT: stats.pearsonr: raise error for complex input
In pearsonr rasie a ValueError for complex data.
Unit test included to show the intended error.

<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
Closes gh-17518

#### What does this implement/fix?
Implement a conditional to check input arrays for complex data.
If the conditional is taken raise ValueError

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
ENH: stats.pearsonr: support complex data
### Describe your issue.

I create a numpy array full of complex numbers, and I compute the correlation of the array with itself. Instead of the expected output of R = 1, some other value is produced.

### Reproducing Code Example

```python
import numpy as np
import scipy.stats

array = np.array([
    -4.21904583e-03-1.53714478e-03j, -4.24663044e-03-1.12832926e-03j,
    -4.26968892e-03-4.87451439e-04j, -6.99917538e-03+3.07376860e-04j,
    0.00000000e+00+0.00000000e+00j, -2.42585590e-02+1.42052459e-02j,
    -5.53404148e-03+4.60188062e-03j, -4.68829482e-03+4.90179019e-03j,
    -7.02331258e-03+8.75908673e-03j, -1.31233383e-01+1.86572484e-01j,
    -4.05137401e-03+6.59972035e-03j, -4.20701822e-03+7.29813816e-03j,
    -3.56487231e-03+6.51759430e-03j, -3.68077200e-03+7.04388575e-03j,
    -8.16459981e-02+1.70084145e-01j, -5.11737898e-03+1.98164995e-02j,
    6.72772914e-04-7.28110367e-05j,  2.13957504e-03-1.82525995e-03j,
    1.60369835e-03-1.54029189e-03j,  8.77788719e-02-8.45568854e-02j,
    1.04277417e-01-9.38854749e-02j,  7.58465696e-03-6.07906563e-03j,
    8.00776452e-03-5.70470615e-03j,  8.36166252e-03-5.14978313e-03j,
    0.00000000e+00+0.00000000e+00j,  0.00000000e+00+0.00000000e+00j,
    0.00000000e+00+0.00000000e+00j,  7.26422461e-03+4.40382166e-04j,
    4.01364547e-03+1.09269127e-03j, -1.99069471e-01-1.20355081e-01j,
    1.56511579e-01+2.59839758e-01j,  9.14046953e-04+5.42262898e-03j,
    -8.37800782e-04+5.67555708e-03j, -3.36561822e-03+7.50108018e-03j,
    -4.22682090e-03+5.36279242e-03j,  5.95438564e-02-3.48209841e-02j,
    -6.77184281e-03+2.10711488e-03j, -4.84293269e-03+3.78698499e-04j,
    -5.13547723e-03-6.86765713e-04j,  4.48392070e-01+1.54568226e-01j,
    -3.17412047e-01-2.35431216e-01j, -2.95731737e-03-3.39078899e-03j,
    -1.95111443e-03-3.77545168e-03j, -2.82719903e-04-1.61393513e-03j,
    7.20241467e-04-1.73515565e-03j, -1.96675563e-01-4.42259734e-02j,
    0.00000000e+00+0.00000000e+00j,  4.84813452e-03+7.60742077e-03j,
    6.31707602e-03+1.51808252e-02j,  2.99277774e-03+1.18667410e-02j,
    5.64640060e-04+1.58372118e-02j, -1.74137347e-03+1.70383706e-02j,
    -5.91398408e-03+2.30008930e-02j, -7.12027831e-03+1.87732435e-02j,
    9.30919156e-02-1.65255887e-01j, -2.09716130e-01+2.30490479e-01j,
    -1.80115101e-02+1.37248240e-02j, -1.85851718e-02+9.23420957e-03j,
    -1.88459965e-02+5.12854226e-03j,  1.09175874e+00-9.17875627e-02j,
    -1.63766142e-02-5.32431671e-03j, -1.24749963e-02-9.63714407e-03j,
    -7.58657222e-03-1.27728267e-02j, -1.99052439e-03-1.35879033e-02j,
    -5.70595470e-01+2.27742231e+00j,  1.24516564e-02-1.21867738e-02j,
    1.82174257e-02-8.67884733e-03j,  2.27204879e-02-3.77097224e-03j,
    2.66143091e-02+2.68683768e-03j,  1.06983372e+00+3.19301893e-01j,
    -6.86033738e-01-4.72910865e-01j,  3.00291320e-02+3.10297521e-02j,
    2.22880055e-02+3.45332319e-02j,  1.61724440e-02+4.04122368e-02j,
    9.78881043e-03+4.96053678e-02j, -6.51085120e-03+5.27227722e-02j,
    -1.76752380e-02+5.26095806e-02j, -3.81856382e-02+6.41735764e-02j,
    0.00000000e+00+0.00000000e+00j, -4.32481463e-02+3.88706950e-02j
])
r, pval = scipy.stats.pearsonr(array, array)
assert np.isclose(r, 1.0), r
```


### Error message

```shell
AssertionError                            Traceback (most recent call last)
Cell In [77], line 47
      4 array = np.array([
      5     -4.21904583e-03-1.53714478e-03j, -4.24663044e-03-1.12832926e-03j,
      6     -4.26968892e-03-4.87451439e-04j, -6.99917538e-03+3.07376860e-04j,
   (...)
     44     0.00000000e+00+0.00000000e+00j, -4.32481463e-02+3.88706950e-02j
     45 ])
     46 r, pval = scipy.stats.pearsonr(array, array)
---> 47 assert np.isclose(r, 1.0), r

AssertionError: (-0.22625405562215772-0.14080706439966875j)
```


### SciPy/NumPy/Python version information

1.9.1 1.23.3 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)

Hints from issue discussion:
I encountered a similar problem in [xarray](https://github.com/pydata/xarray/issues/7340). One of the arrays needs to be conjugated before computing the dot product when computing covariance.
Not much of scipy.stats appears to have been designed with complex data in mind, so I don't think this is a bug. I'll change the title accordingly. 
Shouldn't the function at least raise an exception when passed a complex data type instead of giving the wrong answer?
That would be better, yes. Please submit a PR with that fix and a unit test showing that the intended error (with the intended message) is raised, and I can review and merge it.
@mdhaber Was this addressed in the mentioned pull request or is there still room to submit a PR for this?
This is still waiting on a PR. If you choose to do the complex conjugate approach, please provide a reference for why that is appropriate and demonstrating that the null distribution of the statistic is unchanged (so that p-value and confidence interval are still correct). If we can't show that easily, let's raise a `ValueError` for now. There's probably still time to backport that to version1.10 before release.

Commit message: ENH: stats.pearsonr: issue 17518

In pearsonr rasie a ValueError for complex data.
Unit test included to show the intended error.

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index 2dadadef0c8a..8b8633ecc09b 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -4413,6 +4413,9 @@ def pearsonr(x, y, *, alternative='two-sided'):
     x = np.asarray(x)
     y = np.asarray(y)
 
+    if True in np.iscomplex(x) or True in np.iscomplex(y):
+        raise ValueError('This function does not support complex data')
+
     # If an input is constant, the correlation coefficient is not defined.
     if (x == x[0]).all() or (y == y[0]).all():
         msg = ("An input array is constant; the correlation coefficient "
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 37b93ef35329..4c11a6d5f5fb 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -477,6 +477,11 @@ def test_len1(self):
         y = [2]
         assert_raises(ValueError, stats.pearsonr, x, y)
 
+    def test_complex_data(self):
+        x = [-1j, -2j, -3.0j]
+        y = [-1j, -2j, -3.0j]
+        assert_raises(ValueError, stats.pearsonr, x, y)
+
 
 class TestFisherExact:
     """Some tests to show that fisher_exact() works correctly.


Code review comments:
diff_hunk: @@ -4413,6 +4413,9 @@ def pearsonr(x, y, *, alternative='two-sided'):
     x = np.asarray(x)
     y = np.asarray(y)
 
+    if True in np.iscomplex(x) or True in np.iscomplex(y):
path: scipy/stats/_stats_py.py
text: @user1:
Since we have NumPy arrays, we would usually use `np.any` or the `any` method of the arrays rather than pure-Python.

```suggestion
    if np.iscomplex(x).any() or np.iscomplex(y).any():
```

But that would be slow, because `any` would check every element.

For that reason, I'd suggest checking the dtype:

```suggestion
    if (np.issubdtype(x.dtype, np.complexfloating) 
            or np.issubdtype(y.dtype, np.complexfloating)):
```

diff_hunk: @@ -477,6 +477,11 @@ def test_len1(self):
         y = [2]
         assert_raises(ValueError, stats.pearsonr, x, y)
 
+    def test_complex_data(self):
+        x = [-1j, -2j, -3.0j]
+        y = [-1j, -2j, -3.0j]
+        assert_raises(ValueError, stats.pearsonr, x, y)
path: scipy/stats/tests/test_stats.py
text: @user1:
There are many `ValueErrors` out there. This test should only pass if the one we have in mind is raised.
```suggestion
        message = 'This function does not support complex data'
        with pytest.raises(ValueError, match=message):
            stats.pearsonr(x, y)
```