Merge pull request scipy#13169 from WarrenWeckesser/mstats-pearsonr

BUG: stats: Fix handling of misaligned masks in mstats.pearsonr.

scipy/stats/mstats_basic.py

@@ -355,7 +355,7 @@ def msign(x):
     return ma.filled(np.sign(x), 0)
 
 
-def pearsonr(x,y):
+def pearsonr(x, y):
     """
     Calculates a Pearson correlation coefficient and the p-value for testing
     non-correlation.
@@ -399,24 +399,8 @@ def pearsonr(x,y):
     if df < 0:
         return (masked, masked)
 
-    (mx, my) = (x.mean(), y.mean())
-    (xm, ym) = (x-mx, y-my)
-
-    r_num = ma.add.reduce(xm*ym)
-    r_den = ma.sqrt(ma.dot(xm,xm) * ma.dot(ym,ym))
-    r = r_num / r_den
-    # Presumably, if r > 1, then it is only some small artifact of floating
-    # point arithmetic.
-    r = min(r, 1.0)
-    r = max(r, -1.0)
-
-    if r is masked or abs(r) == 1.0:
-        prob = 0.
-    else:
-        t_squared = (df / ((1.0 - r) * (1.0 + r))) * r * r
-        prob = _betai(0.5*df, 0.5, df/(df + t_squared))
-
-    return r, prob
+    return scipy.stats.stats.pearsonr(ma.masked_array(x, mask=m).compressed(),
+                                      ma.masked_array(y, mask=m).compressed())
 
 
 SpearmanrResult = namedtuple('SpearmanrResult', ('correlation', 'pvalue'))

scipy/stats/tests/test_mstats_basic.py

@@ -201,6 +201,16 @@ def test_pearsonr(self):
         assert_almost_equal(r, np.sqrt(3)/2)
         assert_almost_equal(p, 1.0/3)
 
+    def test_pearsonr_misaligned_mask(self):
+        mx = np.ma.masked_array([1, 2, 3, 4, 5, 6], mask=[0, 1, 0, 0, 0, 0])
+        my = np.ma.masked_array([9, 8, 7, 6, 5, 9], mask=[0, 0, 1, 0, 0, 0])
+        x = np.array([1, 4, 5, 6])
+        y = np.array([9, 6, 5, 9])
+        mr, mp = mstats.pearsonr(mx, my)
+        r, p = stats.pearsonr(x, y)
+        assert_equal(mr, r)
+        assert_equal(mp, p)
+
     def test_spearmanr(self):
         # Tests some computations of Spearman's rho
         (x, y) = ([5.05,6.75,3.21,2.66], [1.65,2.64,2.64,6.95])