Skip to content

Commit

Permalink
BUG: stats: Fix handling of misaligned masks in mstats.pearsonr.
Browse files Browse the repository at this point in the history
After compressing the masked arrays using a common mask for x and
y, the data is passed to stats.pearsonr to do the calculation.

Closes scipygh-3645.
  • Loading branch information
WarrenWeckesser committed Nov 29, 2020
1 parent c3e640c commit b22c8da
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 19 deletions.
22 changes: 3 additions & 19 deletions scipy/stats/mstats_basic.py
Expand Up @@ -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.
Expand Down Expand Up @@ -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'))
Expand Down
10 changes: 10 additions & 0 deletions scipy/stats/tests/test_mstats_basic.py
Expand Up @@ -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])
Expand Down

0 comments on commit b22c8da

Please sign in to comment.