master
Per.Andreas.Brodtkorb 11 years ago
parent b9031f8a22
commit c070dd4dc4

@ -14,11 +14,13 @@ import wafo.stats.mstats as mstats
from wafo import stats from wafo import stats
from numpy.testing import TestCase, run_module_suite from numpy.testing import TestCase, run_module_suite
from numpy.ma.testutils import (assert_equal, assert_almost_equal, from numpy.ma.testutils import (assert_equal, assert_almost_equal,
assert_array_almost_equal, assert_array_almost_equal_nulp, assert_, assert_array_almost_equal,
assert_array_almost_equal_nulp, assert_,
assert_allclose, assert_raises) assert_allclose, assert_raises)
class TestMquantiles(TestCase): class TestMquantiles(TestCase):
def test_mquantiles_limit_keyword(self): def test_mquantiles_limit_keyword(self):
# Regression test for Trac ticket #867 # Regression test for Trac ticket #867
data = np.array([[6., 7., 1.], data = np.array([[6., 7., 1.],
@ -40,68 +42,70 @@ class TestMquantiles(TestCase):
class TestGMean(TestCase): class TestGMean(TestCase):
def test_1D(self): def test_1D(self):
a = (1,2,3,4) a = (1, 2, 3, 4)
actual = mstats.gmean(a) actual = mstats.gmean(a)
desired = np.power(1*2*3*4,1./4.) desired = np.power(1 * 2 * 3 * 4, 1. / 4.)
assert_almost_equal(actual, desired,decimal=14) assert_almost_equal(actual, desired, decimal=14)
desired1 = mstats.gmean(a,axis=-1) desired1 = mstats.gmean(a, axis=-1)
assert_almost_equal(actual, desired1, decimal=14) assert_almost_equal(actual, desired1, decimal=14)
assert_(not isinstance(desired1, ma.MaskedArray)) assert_(not isinstance(desired1, ma.MaskedArray))
a = ma.array((1,2,3,4),mask=(0,0,0,1)) a = ma.array((1, 2, 3, 4), mask=(0, 0, 0, 1))
actual = mstats.gmean(a) actual = mstats.gmean(a)
desired = np.power(1*2*3,1./3.) desired = np.power(1 * 2 * 3, 1. / 3.)
assert_almost_equal(actual, desired,decimal=14) assert_almost_equal(actual, desired, decimal=14)
desired1 = mstats.gmean(a,axis=-1) desired1 = mstats.gmean(a, axis=-1)
assert_almost_equal(actual, desired1, decimal=14) assert_almost_equal(actual, desired1, decimal=14)
def test_2D(self): def test_2D(self):
a = ma.array(((1,2,3,4),(1,2,3,4),(1,2,3,4)), a = ma.array(((1, 2, 3, 4), (1, 2, 3, 4), (1, 2, 3, 4)),
mask=((0,0,0,0),(1,0,0,1),(0,1,1,0))) mask=((0, 0, 0, 0), (1, 0, 0, 1), (0, 1, 1, 0)))
actual = mstats.gmean(a) actual = mstats.gmean(a)
desired = np.array((1,2,3,4)) desired = np.array((1, 2, 3, 4))
assert_array_almost_equal(actual, desired, decimal=14) assert_array_almost_equal(actual, desired, decimal=14)
desired1 = mstats.gmean(a,axis=0) desired1 = mstats.gmean(a, axis=0)
assert_array_almost_equal(actual, desired1, decimal=14) assert_array_almost_equal(actual, desired1, decimal=14)
actual = mstats.gmean(a, -1) actual = mstats.gmean(a, -1)
desired = ma.array((np.power(1*2*3*4,1./4.), desired = ma.array((np.power(1 * 2 * 3 * 4, 1. / 4.),
np.power(2*3,1./2.), np.power(2 * 3, 1. / 2.),
np.power(1*4,1./2.))) np.power(1 * 4, 1. / 2.)))
assert_array_almost_equal(actual, desired, decimal=14) assert_array_almost_equal(actual, desired, decimal=14)
class TestHMean(TestCase): class TestHMean(TestCase):
def test_1D(self): def test_1D(self):
a = (1,2,3,4) a = (1, 2, 3, 4)
actual = mstats.hmean(a) actual = mstats.hmean(a)
desired = 4. / (1./1 + 1./2 + 1./3 + 1./4) desired = 4. / (1. / 1 + 1. / 2 + 1. / 3 + 1. / 4)
assert_almost_equal(actual, desired, decimal=14) assert_almost_equal(actual, desired, decimal=14)
desired1 = mstats.hmean(ma.array(a),axis=-1) desired1 = mstats.hmean(ma.array(a), axis=-1)
assert_almost_equal(actual, desired1, decimal=14) assert_almost_equal(actual, desired1, decimal=14)
a = ma.array((1,2,3,4),mask=(0,0,0,1)) a = ma.array((1, 2, 3, 4), mask=(0, 0, 0, 1))
actual = mstats.hmean(a) actual = mstats.hmean(a)
desired = 3. / (1./1 + 1./2 + 1./3) desired = 3. / (1. / 1 + 1. / 2 + 1. / 3)
assert_almost_equal(actual, desired,decimal=14) assert_almost_equal(actual, desired, decimal=14)
desired1 = mstats.hmean(a,axis=-1) desired1 = mstats.hmean(a, axis=-1)
assert_almost_equal(actual, desired1, decimal=14) assert_almost_equal(actual, desired1, decimal=14)
def test_2D(self): def test_2D(self):
a = ma.array(((1,2,3,4),(1,2,3,4),(1,2,3,4)), a = ma.array(((1, 2, 3, 4), (1, 2, 3, 4), (1, 2, 3, 4)),
mask=((0,0,0,0),(1,0,0,1),(0,1,1,0))) mask=((0, 0, 0, 0), (1, 0, 0, 1), (0, 1, 1, 0)))
actual = mstats.hmean(a) actual = mstats.hmean(a)
desired = ma.array((1,2,3,4)) desired = ma.array((1, 2, 3, 4))
assert_array_almost_equal(actual, desired, decimal=14) assert_array_almost_equal(actual, desired, decimal=14)
actual1 = mstats.hmean(a,axis=-1) actual1 = mstats.hmean(a, axis=-1)
desired = (4./(1/1.+1/2.+1/3.+1/4.), desired = (4. / (1 / 1. + 1 / 2. + 1 / 3. + 1 / 4.),
2./(1/2.+1/3.), 2. / (1 / 2. + 1 / 3.),
2./(1/1.+1/4.) 2. / (1 / 1. + 1 / 4.)
) )
assert_array_almost_equal(actual1, desired, decimal=14) assert_array_almost_equal(actual1, desired, decimal=14)
@ -112,18 +116,24 @@ class TestRanking(TestCase):
TestCase.__init__(self, *args, **kwargs) TestCase.__init__(self, *args, **kwargs)
def test_ranking(self): def test_ranking(self):
x = ma.array([0,1,1,1,2,3,4,5,5,6,]) x = ma.array([0, 1, 1, 1, 2, 3, 4, 5, 5, 6, ])
assert_almost_equal(mstats.rankdata(x),[1,3,3,3,5,6,7,8.5,8.5,10]) assert_almost_equal(
x[[3,4]] = masked mstats.rankdata(x), [1, 3, 3, 3, 5, 6, 7, 8.5, 8.5, 10])
assert_almost_equal(mstats.rankdata(x),[1,2.5,2.5,0,0,4,5,6.5,6.5,8]) x[[3, 4]] = masked
assert_almost_equal(mstats.rankdata(x,use_missing=True), assert_almost_equal(
[1,2.5,2.5,4.5,4.5,4,5,6.5,6.5,8]) mstats.rankdata(x), [1, 2.5, 2.5, 0, 0, 4, 5, 6.5, 6.5, 8])
x = ma.array([0,1,5,1,2,4,3,5,1,6,]) assert_almost_equal(mstats.rankdata(x, use_missing=True),
assert_almost_equal(mstats.rankdata(x),[1,3,8.5,3,5,7,6,8.5,3,10]) [1, 2.5, 2.5, 4.5, 4.5, 4, 5, 6.5, 6.5, 8])
x = ma.array([[0,1,1,1,2], [3,4,5,5,6,]]) x = ma.array([0, 1, 5, 1, 2, 4, 3, 5, 1, 6, ])
assert_almost_equal(mstats.rankdata(x),[[1,3,3,3,5],[6,7,8.5,8.5,10]]) assert_almost_equal(
assert_almost_equal(mstats.rankdata(x,axis=1),[[1,3,3,3,5],[1,2,3.5,3.5,5]]) mstats.rankdata(x), [1, 3, 8.5, 3, 5, 7, 6, 8.5, 3, 10])
assert_almost_equal(mstats.rankdata(x,axis=0),[[1,1,1,1,1],[2,2,2,2,2,]]) x = ma.array([[0, 1, 1, 1, 2], [3, 4, 5, 5, 6, ]])
assert_almost_equal(
mstats.rankdata(x), [[1, 3, 3, 3, 5], [6, 7, 8.5, 8.5, 10]])
assert_almost_equal(
mstats.rankdata(x, axis=1), [[1, 3, 3, 3, 5], [1, 2, 3.5, 3.5, 5]])
assert_almost_equal(
mstats.rankdata(x, axis=0), [[1, 1, 1, 1, 1], [2, 2, 2, 2, 2, ]])
class TestCorr(TestCase): class TestCorr(TestCase):
@ -148,72 +158,73 @@ class TestCorr(TestCase):
x1 = ma.array([-1.0, 0.0, 1.0]) x1 = ma.array([-1.0, 0.0, 1.0])
y1 = ma.array([0, 0, 3]) y1 = ma.array([0, 0, 3])
r, p = mstats.pearsonr(x1, y1) r, p = mstats.pearsonr(x1, y1)
assert_almost_equal(r, np.sqrt(3)/2) assert_almost_equal(r, np.sqrt(3) / 2)
assert_almost_equal(p, 1.0/3) assert_almost_equal(p, 1.0 / 3)
# (x2, y2) have the same unmasked data as (x1, y1). # (x2, y2) have the same unmasked data as (x1, y1).
mask = [False, False, False, True] mask = [False, False, False, True]
x2 = ma.array([-1.0, 0.0, 1.0, 99.0], mask=mask) x2 = ma.array([-1.0, 0.0, 1.0, 99.0], mask=mask)
y2 = ma.array([0, 0, 3, -1], mask=mask) y2 = ma.array([0, 0, 3, -1], mask=mask)
r, p = mstats.pearsonr(x2, y2) r, p = mstats.pearsonr(x2, y2)
assert_almost_equal(r, np.sqrt(3)/2) assert_almost_equal(r, np.sqrt(3) / 2)
assert_almost_equal(p, 1.0/3) assert_almost_equal(p, 1.0 / 3)
def test_spearmanr(self): def test_spearmanr(self):
# Tests some computations of Spearman's rho # 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]) (x, y) = ([5.05, 6.75, 3.21, 2.66], [1.65, 2.64, 2.64, 6.95])
assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555) assert_almost_equal(mstats.spearmanr(x, y)[0], -0.6324555)
(x, y) = ([5.05,6.75,3.21,2.66,np.nan],[1.65,2.64,2.64,6.95,np.nan]) (x, y) = ([5.05, 6.75, 3.21, 2.66, np.nan],
[1.65, 2.64, 2.64, 6.95, np.nan])
(x, y) = (ma.fix_invalid(x), ma.fix_invalid(y)) (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555) assert_almost_equal(mstats.spearmanr(x, y)[0], -0.6324555)
x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1,
1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7] 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7]
y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6,
0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4] 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4]
assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299) assert_almost_equal(mstats.spearmanr(x, y)[0], 0.6887299)
x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, x = [2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1,
1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7, np.nan] 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7, np.nan]
y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, y = [22.6, 8.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6,
0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4, np.nan] 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4, np.nan]
(x, y) = (ma.fix_invalid(x), ma.fix_invalid(y)) (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y))
assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299) assert_almost_equal(mstats.spearmanr(x, y)[0], 0.6887299)
def test_kendalltau(self): def test_kendalltau(self):
# Tests some computations of Kendall's tau # Tests some computations of Kendall's tau
x = ma.fix_invalid([5.05, 6.75, 3.21, 2.66,np.nan]) x = ma.fix_invalid([5.05, 6.75, 3.21, 2.66, np.nan])
y = ma.fix_invalid([1.65, 26.5, -5.93, 7.96, np.nan]) y = ma.fix_invalid([1.65, 26.5, -5.93, 7.96, np.nan])
z = ma.fix_invalid([1.65, 2.64, 2.64, 6.95, np.nan]) z = ma.fix_invalid([1.65, 2.64, 2.64, 6.95, np.nan])
assert_almost_equal(np.asarray(mstats.kendalltau(x,y)), assert_almost_equal(np.asarray(mstats.kendalltau(x, y)),
[+0.3333333,0.4969059]) [+0.3333333, 0.4969059])
assert_almost_equal(np.asarray(mstats.kendalltau(x,z)), assert_almost_equal(np.asarray(mstats.kendalltau(x, z)),
[-0.5477226,0.2785987]) [-0.5477226, 0.2785987])
# #
x = ma.fix_invalid([0, 0, 0, 0,20,20, 0,60, 0,20, x = ma.fix_invalid([0, 0, 0, 0, 20, 20, 0, 60, 0, 20,
10,10, 0,40, 0,20, 0, 0, 0, 0, 0, np.nan]) 10, 10, 0, 40, 0, 20, 0, 0, 0, 0, 0, np.nan])
y = ma.fix_invalid([0,80,80,80,10,33,60, 0,67,27, y = ma.fix_invalid([0, 80, 80, 80, 10, 33, 60, 0, 67, 27,
25,80,80,80,80,80,80, 0,10,45, np.nan, 0]) 25, 80, 80, 80, 80, 80, 80, 0, 10, 45, np.nan, 0])
result = mstats.kendalltau(x,y) result = mstats.kendalltau(x, y)
assert_almost_equal(np.asarray(result), [-0.1585188, 0.4128009]) assert_almost_equal(np.asarray(result), [-0.1585188, 0.4128009])
def test_kendalltau_seasonal(self): def test_kendalltau_seasonal(self):
# Tests the seasonal Kendall tau. # Tests the seasonal Kendall tau.
x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], x = [[nan, nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
[4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
[3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], [3, 2, 5, 6, 18, 4, 9, 1, 1, nan, 1, 1, nan],
[nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] [nan, 6, 11, 4, 17, nan, 6, 1, 1, 2, 5, 1, 1]]
x = ma.fix_invalid(x).T x = ma.fix_invalid(x).T
output = mstats.kendalltau_seasonal(x) output = mstats.kendalltau_seasonal(x)
assert_almost_equal(output['global p-value (indep)'], 0.008, 3) assert_almost_equal(output['global p-value (indep)'], 0.008, 3)
assert_almost_equal(output['seasonal p-value'].round(2), assert_almost_equal(output['seasonal p-value'].round(2),
[0.18,0.53,0.20,0.04]) [0.18, 0.53, 0.20, 0.04])
def test_pointbiserial(self): def test_pointbiserial(self):
x = [1,0,1,1,1,1,0,1,0,0,0,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0, x = [1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
0,0,0,0,1,-1] 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, -1]
y = [14.8,13.8,12.4,10.1,7.1,6.1,5.8,4.6,4.3,3.5,3.3,3.2,3.0, y = [14.8, 13.8, 12.4, 10.1, 7.1, 6.1, 5.8, 4.6, 4.3, 3.5, 3.3, 3.2, 3.0,
2.8,2.8,2.5,2.4,2.3,2.1,1.7,1.7,1.5,1.3,1.3,1.2,1.2,1.1, 2.8, 2.8, 2.5, 2.4, 2.3, 2.1, 1.7, 1.7, 1.5, 1.3, 1.3, 1.2, 1.2, 1.1,
0.8,0.7,0.6,0.5,0.2,0.2,0.1,np.nan] 0.8, 0.7, 0.6, 0.5, 0.2, 0.2, 0.1, np.nan]
assert_almost_equal(mstats.pointbiserialr(x, y)[0], 0.36149, 5) assert_almost_equal(mstats.pointbiserialr(x, y)[0], 0.36149, 5)
@ -221,68 +232,70 @@ class TestTrimming(TestCase):
def test_trim(self): def test_trim(self):
a = ma.arange(10) a = ma.arange(10)
assert_equal(mstats.trim(a), [0,1,2,3,4,5,6,7,8,9]) assert_equal(mstats.trim(a), [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
a = ma.arange(10) a = ma.arange(10)
assert_equal(mstats.trim(a,(2,8)), [None,None,2,3,4,5,6,7,8,None]) assert_equal(
mstats.trim(a, (2, 8)), [None, None, 2, 3, 4, 5, 6, 7, 8, None])
a = ma.arange(10) a = ma.arange(10)
assert_equal(mstats.trim(a,limits=(2,8),inclusive=(False,False)), assert_equal(mstats.trim(a, limits=(2, 8), inclusive=(False, False)),
[None,None,None,3,4,5,6,7,None,None]) [None, None, None, 3, 4, 5, 6, 7, None, None])
a = ma.arange(10) a = ma.arange(10)
assert_equal(mstats.trim(a,limits=(0.1,0.2),relative=True), assert_equal(mstats.trim(a, limits=(0.1, 0.2), relative=True),
[None,1,2,3,4,5,6,7,None,None]) [None, 1, 2, 3, 4, 5, 6, 7, None, None])
a = ma.arange(12) a = ma.arange(12)
a[[0,-1]] = a[5] = masked a[[0, -1]] = a[5] = masked
assert_equal(mstats.trim(a,(2,8)), assert_equal(mstats.trim(a, (2, 8)),
[None,None,2,3,4,None,6,7,8,None,None,None]) [None, None, 2, 3, 4, None, 6, 7, 8, None, None, None])
x = ma.arange(100).reshape(10,10) x = ma.arange(100).reshape(10, 10)
trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=None) trimx = mstats.trim(x, (0.1, 0.2), relative=True, axis=None)
assert_equal(trimx._mask.ravel(),[1]*10+[0]*70+[1]*20) assert_equal(trimx._mask.ravel(), [1] * 10 + [0] * 70 + [1] * 20)
trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=0) trimx = mstats.trim(x, (0.1, 0.2), relative=True, axis=0)
assert_equal(trimx._mask.ravel(),[1]*10+[0]*70+[1]*20) assert_equal(trimx._mask.ravel(), [1] * 10 + [0] * 70 + [1] * 20)
trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=-1) trimx = mstats.trim(x, (0.1, 0.2), relative=True, axis=-1)
assert_equal(trimx._mask.T.ravel(),[1]*10+[0]*70+[1]*20) assert_equal(trimx._mask.T.ravel(), [1] * 10 + [0] * 70 + [1] * 20)
x = ma.arange(110).reshape(11,10) x = ma.arange(110).reshape(11, 10)
x[1] = masked x[1] = masked
trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=None) trimx = mstats.trim(x, (0.1, 0.2), relative=True, axis=None)
assert_equal(trimx._mask.ravel(),[1]*20+[0]*70+[1]*20) assert_equal(trimx._mask.ravel(), [1] * 20 + [0] * 70 + [1] * 20)
trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=0) trimx = mstats.trim(x, (0.1, 0.2), relative=True, axis=0)
assert_equal(trimx._mask.ravel(),[1]*20+[0]*70+[1]*20) assert_equal(trimx._mask.ravel(), [1] * 20 + [0] * 70 + [1] * 20)
trimx = mstats.trim(x.T,(0.1,0.2),relative=True,axis=-1) trimx = mstats.trim(x.T, (0.1, 0.2), relative=True, axis=-1)
assert_equal(trimx.T._mask.ravel(),[1]*20+[0]*70+[1]*20) assert_equal(trimx.T._mask.ravel(), [1] * 20 + [0] * 70 + [1] * 20)
def test_trim_old(self): def test_trim_old(self):
x = ma.arange(100) x = ma.arange(100)
assert_equal(mstats.trimboth(x).count(), 60) assert_equal(mstats.trimboth(x).count(), 60)
assert_equal(mstats.trimtail(x,tail='r').count(), 80) assert_equal(mstats.trimtail(x, tail='r').count(), 80)
x[50:70] = masked x[50:70] = masked
trimx = mstats.trimboth(x) trimx = mstats.trimboth(x)
assert_equal(trimx.count(), 48) assert_equal(trimx.count(), 48)
assert_equal(trimx._mask, [1]*16 + [0]*34 + [1]*20 + [0]*14 + [1]*16) assert_equal(
trimx._mask, [1] * 16 + [0] * 34 + [1] * 20 + [0] * 14 + [1] * 16)
x._mask = nomask x._mask = nomask
x.shape = (10,10) x.shape = (10, 10)
assert_equal(mstats.trimboth(x).count(), 60) assert_equal(mstats.trimboth(x).count(), 60)
assert_equal(mstats.trimtail(x).count(), 80) assert_equal(mstats.trimtail(x).count(), 80)
def test_trimmedmean(self): def test_trimmedmean(self):
data = ma.array([77, 87, 88,114,151,210,219,246,253,262, data = ma.array([77, 87, 88, 114, 151, 210, 219, 246, 253, 262,
296,299,306,376,428,515,666,1310,2611]) 296, 299, 306, 376, 428, 515, 666, 1310, 2611])
assert_almost_equal(mstats.trimmed_mean(data,0.1), 343, 0) assert_almost_equal(mstats.trimmed_mean(data, 0.1), 343, 0)
assert_almost_equal(mstats.trimmed_mean(data,(0.1,0.1)), 343, 0) assert_almost_equal(mstats.trimmed_mean(data, (0.1, 0.1)), 343, 0)
assert_almost_equal(mstats.trimmed_mean(data,(0.2,0.2)), 283, 0) assert_almost_equal(mstats.trimmed_mean(data, (0.2, 0.2)), 283, 0)
def test_trimmed_stde(self): def test_trimmed_stde(self):
data = ma.array([77, 87, 88,114,151,210,219,246,253,262, data = ma.array([77, 87, 88, 114, 151, 210, 219, 246, 253, 262,
296,299,306,376,428,515,666,1310,2611]) 296, 299, 306, 376, 428, 515, 666, 1310, 2611])
assert_almost_equal(mstats.trimmed_stde(data,(0.2,0.2)), 56.13193, 5) assert_almost_equal(mstats.trimmed_stde(data, (0.2, 0.2)), 56.13193, 5)
assert_almost_equal(mstats.trimmed_stde(data,0.2), 56.13193, 5) assert_almost_equal(mstats.trimmed_stde(data, 0.2), 56.13193, 5)
def test_winsorization(self): def test_winsorization(self):
data = ma.array([77, 87, 88,114,151,210,219,246,253,262, data = ma.array([77, 87, 88, 114, 151, 210, 219, 246, 253, 262,
296,299,306,376,428,515,666,1310,2611]) 296, 299, 306, 376, 428, 515, 666, 1310, 2611])
assert_almost_equal(mstats.winsorize(data,(0.2,0.2)).var(ddof=1), assert_almost_equal(mstats.winsorize(data, (0.2, 0.2)).var(ddof=1),
21551.4, 1) 21551.4, 1)
data[5] = masked data[5] = masked
winsorized = mstats.winsorize(data) winsorized = mstats.winsorize(data)
@ -297,7 +310,7 @@ class TestMoments(TestCase):
# http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml # http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml
# http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml # http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml
# Note that both test cases came from here. # Note that both test cases came from here.
testcase = [1,2,3,4] testcase = [1, 2, 3, 4]
testmathworks = ma.fix_invalid([1.165, 0.6268, 0.0751, 0.3516, -0.6965, testmathworks = ma.fix_invalid([1.165, 0.6268, 0.0751, 0.3516, -0.6965,
np.nan]) np.nan])
testcase_2d = ma.array( testcase_2d = ma.array(
@ -313,44 +326,44 @@ class TestMoments(TestCase):
[False, False, True, False, False]], dtype=np.bool)) [False, False, True, False, False]], dtype=np.bool))
def test_moment(self): def test_moment(self):
y = mstats.moment(self.testcase,1) y = mstats.moment(self.testcase, 1)
assert_almost_equal(y,0.0,10) assert_almost_equal(y, 0.0, 10)
y = mstats.moment(self.testcase,2) y = mstats.moment(self.testcase, 2)
assert_almost_equal(y,1.25) assert_almost_equal(y, 1.25)
y = mstats.moment(self.testcase,3) y = mstats.moment(self.testcase, 3)
assert_almost_equal(y,0.0) assert_almost_equal(y, 0.0)
y = mstats.moment(self.testcase,4) y = mstats.moment(self.testcase, 4)
assert_almost_equal(y,2.5625) assert_almost_equal(y, 2.5625)
def test_variation(self): def test_variation(self):
y = mstats.variation(self.testcase) y = mstats.variation(self.testcase)
assert_almost_equal(y,0.44721359549996, 10) assert_almost_equal(y, 0.44721359549996, 10)
def test_skewness(self): def test_skewness(self):
y = mstats.skew(self.testmathworks) y = mstats.skew(self.testmathworks)
assert_almost_equal(y,-0.29322304336607,10) assert_almost_equal(y, -0.29322304336607, 10)
y = mstats.skew(self.testmathworks,bias=0) y = mstats.skew(self.testmathworks, bias=0)
assert_almost_equal(y,-0.437111105023940,10) assert_almost_equal(y, -0.437111105023940, 10)
y = mstats.skew(self.testcase) y = mstats.skew(self.testcase)
assert_almost_equal(y,0.0,10) assert_almost_equal(y, 0.0, 10)
def test_kurtosis(self): def test_kurtosis(self):
# Set flags for axis = 0 and fisher=0 (Pearson's definition of kurtosis # Set flags for axis = 0 and fisher=0 (Pearson's definition of kurtosis
# for compatibility with Matlab) # for compatibility with Matlab)
y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1) y = mstats.kurtosis(self.testmathworks, 0, fisher=0, bias=1)
assert_almost_equal(y, 2.1658856802973,10) assert_almost_equal(y, 2.1658856802973, 10)
# Note that MATLAB has confusing docs for the following case # Note that MATLAB has confusing docs for the following case
# kurtosis(x,0) gives an unbiased estimate of Pearson's skewness # kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
# kurtosis(x) gives a biased estimate of Fisher's skewness (Pearson-3) # kurtosis(x) gives a biased estimate of Fisher's skewness (Pearson-3)
# The MATLAB docs imply that both should give Fisher's # The MATLAB docs imply that both should give Fisher's
y = mstats.kurtosis(self.testmathworks,fisher=0, bias=0) y = mstats.kurtosis(self.testmathworks, fisher=0, bias=0)
assert_almost_equal(y, 3.663542721189047,10) assert_almost_equal(y, 3.663542721189047, 10)
y = mstats.kurtosis(self.testcase,0,0) y = mstats.kurtosis(self.testcase, 0, 0)
assert_almost_equal(y,1.64) assert_almost_equal(y, 1.64)
# test that kurtosis works on multidimensional masked arrays # test that kurtosis works on multidimensional masked arrays
correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0., correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385, 0.,
-1.26979517952]), - 1.26979517952]),
mask=np.array([False, False, False, True, mask=np.array([False, False, False, True,
False], dtype=np.bool)) False], dtype=np.bool))
assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1), assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
@ -373,34 +386,40 @@ class TestMoments(TestCase):
stats.kurtosis(self.testcase_2d[2, :])) stats.kurtosis(self.testcase_2d[2, :]))
def test_mode(self): def test_mode(self):
a1 = [0,0,0,1,1,1,2,3,3,3,3,4,5,6,7] a1 = [0, 0, 0, 1, 1, 1, 2, 3, 3, 3, 3, 4, 5, 6, 7]
a2 = np.reshape(a1, (3,5)) a2 = np.reshape(a1, (3, 5))
a3 = np.array([1,2,3,4,5,6]) a3 = np.array([1, 2, 3, 4, 5, 6])
a4 = np.reshape(a3, (3,2)) a4 = np.reshape(a3, (3, 2))
ma1 = ma.masked_where(ma.array(a1) > 2, a1) ma1 = ma.masked_where(ma.array(a1) > 2, a1)
ma2 = ma.masked_where(a2 > 2, a2) ma2 = ma.masked_where(a2 > 2, a2)
ma3 = ma.masked_where(a3 < 2, a3) ma3 = ma.masked_where(a3 < 2, a3)
ma4 = ma.masked_where(ma.array(a4) < 2, a4) ma4 = ma.masked_where(ma.array(a4) < 2, a4)
assert_equal(mstats.mode(a1, axis=None), (3,4)) assert_equal(mstats.mode(a1, axis=None), (3, 4))
assert_equal(mstats.mode(a1, axis=0), (3,4)) assert_equal(mstats.mode(a1, axis=0), (3, 4))
assert_equal(mstats.mode(ma1, axis=None), (0,3)) assert_equal(mstats.mode(ma1, axis=None), (0, 3))
assert_equal(mstats.mode(a2, axis=None), (3,4)) assert_equal(mstats.mode(a2, axis=None), (3, 4))
assert_equal(mstats.mode(ma2, axis=None), (0,3)) assert_equal(mstats.mode(ma2, axis=None), (0, 3))
assert_equal(mstats.mode(a3, axis=None), (1,1)) assert_equal(mstats.mode(a3, axis=None), (1, 1))
assert_equal(mstats.mode(ma3, axis=None), (2,1)) assert_equal(mstats.mode(ma3, axis=None), (2, 1))
assert_equal(mstats.mode(a2, axis=0), ([[0,0,0,1,1]], [[1,1,1,1,1]])) assert_equal(
assert_equal(mstats.mode(ma2, axis=0), ([[0,0,0,1,1]], [[1,1,1,1,1]])) mstats.mode(a2, axis=0), ([[0, 0, 0, 1, 1]], [[1, 1, 1, 1, 1]]))
assert_equal(mstats.mode(a2, axis=-1), ([[0],[3],[3]], [[3],[3],[1]])) assert_equal(
assert_equal(mstats.mode(ma2, axis=-1), ([[0],[1],[0]], [[3],[1],[0]])) mstats.mode(ma2, axis=0), ([[0, 0, 0, 1, 1]], [[1, 1, 1, 1, 1]]))
assert_equal(mstats.mode(ma4, axis=0), ([[3,2]], [[1,1]])) assert_equal(
assert_equal(mstats.mode(ma4, axis=-1), ([[2],[3],[5]], [[1],[1],[1]])) mstats.mode(a2, axis=-1), ([[0], [3], [3]], [[3], [3], [1]]))
assert_equal(
mstats.mode(ma2, axis=-1), ([[0], [1], [0]], [[3], [1], [0]]))
assert_equal(mstats.mode(ma4, axis=0), ([[3, 2]], [[1, 1]]))
assert_equal(
mstats.mode(ma4, axis=-1), ([[2], [3], [5]], [[1], [1], [1]]))
class TestPercentile(TestCase): class TestPercentile(TestCase):
def setUp(self): def setUp(self):
self.a1 = [3,4,5,10,-3,-5,6] self.a1 = [3, 4, 5, 10, -3, -5, 6]
self.a2 = [3,-6,-2,8,7,4,2,1] self.a2 = [3, -6, -2, 8, 7, 4, 2, 1]
self.a3 = [3.,4,5,10,-3,-5,-6,7.0] self.a3 = [3., 4, 5, 10, -3, -5, -6, 7.0]
def test_percentile(self): def test_percentile(self):
x = np.arange(8) * 0.5 x = np.arange(8) * 0.5
@ -414,27 +433,28 @@ class TestPercentile(TestCase):
[4, 4, 3], [4, 4, 3],
[1, 1, 1], [1, 1, 1],
[1, 1, 1]]) [1, 1, 1]])
assert_equal(mstats.scoreatpercentile(x,50), [1,1,1]) assert_equal(mstats.scoreatpercentile(x, 50), [1, 1, 1])
class TestVariability(TestCase): class TestVariability(TestCase):
""" Comparison numbers are found using R v.1.5.1 """ Comparison numbers are found using R v.1.5.1
note that length(testcase) = 4 note that length(testcase) = 4
""" """
testcase = ma.fix_invalid([1,2,3,4,np.nan]) testcase = ma.fix_invalid([1, 2, 3, 4, np.nan])
def test_signaltonoise(self): def test_signaltonoise(self):
# This is not in R, so used: # This is not in R, so used:
# mean(testcase, axis=0) / (sqrt(var(testcase)*3/4)) # mean(testcase, axis=0) / (sqrt(var(testcase)*3/4))
y = mstats.signaltonoise(self.testcase) y = mstats.signaltonoise(self.testcase)
assert_almost_equal(y,2.236067977) assert_almost_equal(y, 2.236067977)
def test_sem(self): def test_sem(self):
# This is not in R, so used: sqrt(var(testcase)*3/4) / sqrt(3) # This is not in R, so used: sqrt(var(testcase)*3/4) / sqrt(3)
y = mstats.sem(self.testcase) y = mstats.sem(self.testcase)
assert_almost_equal(y, 0.6454972244) assert_almost_equal(y, 0.6454972244)
n = self.testcase.count() n = self.testcase.count()
assert_allclose(mstats.sem(self.testcase, ddof=0) * np.sqrt(n/(n-2)), assert_allclose(mstats.sem(self.testcase, ddof=0) * np.sqrt(n / (n - 2)),
mstats.sem(self.testcase, ddof=2)) mstats.sem(self.testcase, ddof=2))
def test_zmap(self): def test_zmap(self):
@ -458,41 +478,41 @@ class TestVariability(TestCase):
class TestMisc(TestCase): class TestMisc(TestCase):
def test_obrientransform(self): def test_obrientransform(self):
args = [[5]*5+[6]*11+[7]*9+[8]*3+[9]*2+[10]*2, args = [[5] * 5 + [6] * 11 + [7] * 9 + [8] * 3 + [9] * 2 + [10] * 2,
[6]+[7]*2+[8]*4+[9]*9+[10]*16] [6] + [7] * 2 + [8] * 4 + [9] * 9 + [10] * 16]
result = [5*[3.1828]+11*[0.5591]+9*[0.0344]+3*[1.6086]+2*[5.2817]+2*[11.0538], result = [5 * [3.1828] + 11 * [0.5591] + 9 * [0.0344] + 3 * [1.6086] + 2 * [5.2817] + 2 * [11.0538],
[10.4352]+2*[4.8599]+4*[1.3836]+9*[0.0061]+16*[0.7277]] [10.4352] + 2 * [4.8599] + 4 * [1.3836] + 9 * [0.0061] + 16 * [0.7277]]
assert_almost_equal(np.round(mstats.obrientransform(*args).T,4), assert_almost_equal(np.round(mstats.obrientransform(*args).T, 4),
result,4) result, 4)
def test_kstwosamp(self): def test_kstwosamp(self):
x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], x = [[nan, nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
[4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
[3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], [3, 2, 5, 6, 18, 4, 9, 1, 1, nan, 1, 1, nan],
[nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] [nan, 6, 11, 4, 17, nan, 6, 1, 1, 2, 5, 1, 1]]
x = ma.fix_invalid(x).T x = ma.fix_invalid(x).T
(winter,spring,summer,fall) = x.T (winter, spring, summer, fall) = x.T
assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring),4), assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring), 4),
(0.1818,0.9892)) (0.1818, 0.9892))
assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'g'),4), assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring, 'g'), 4),
(0.1469,0.7734)) (0.1469, 0.7734))
assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'l'),4), assert_almost_equal(np.round(mstats.ks_twosamp(winter, spring, 'l'), 4),
(0.1818,0.6744)) (0.1818, 0.6744))
def test_friedmanchisq(self): def test_friedmanchisq(self):
# No missing values # No missing values
args = ([9.0,9.5,5.0,7.5,9.5,7.5,8.0,7.0,8.5,6.0], args = ([9.0, 9.5, 5.0, 7.5, 9.5, 7.5, 8.0, 7.0, 8.5, 6.0],
[7.0,6.5,7.0,7.5,5.0,8.0,6.0,6.5,7.0,7.0], [7.0, 6.5, 7.0, 7.5, 5.0, 8.0, 6.0, 6.5, 7.0, 7.0],
[6.0,8.0,4.0,6.0,7.0,6.5,6.0,4.0,6.5,3.0]) [6.0, 8.0, 4.0, 6.0, 7.0, 6.5, 6.0, 4.0, 6.5, 3.0])
result = mstats.friedmanchisquare(*args) result = mstats.friedmanchisquare(*args)
assert_almost_equal(result[0], 10.4737, 4) assert_almost_equal(result[0], 10.4737, 4)
assert_almost_equal(result[1], 0.005317, 6) assert_almost_equal(result[1], 0.005317, 6)
# Missing values # Missing values
x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], x = [[nan, nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1],
[4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], [4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3],
[3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], [3, 2, 5, 6, 18, 4, 9, 1, 1, nan, 1, 1, nan],
[nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] [nan, 6, 11, 4, 17, nan, 6, 1, 1, 2, 5, 1, 1]]
x = ma.fix_invalid(x) x = ma.fix_invalid(x)
result = mstats.friedmanchisquare(*x) result = mstats.friedmanchisquare(*x)
assert_almost_equal(result[0], 2.0156, 4) assert_almost_equal(result[0], 2.0156, 4)
@ -519,7 +539,7 @@ def test_plotting_positions():
class TestNormalitytests(): class TestNormalitytests():
def test_vs_nonmasked(self): def test_vs_nonmasked(self):
x = np.array((-2,-1,0,1,2,3)*4)**2 x = np.array((-2, -1, 0, 1, 2, 3) * 4) ** 2
assert_array_almost_equal(mstats.normaltest(x), stats.normaltest(x)) assert_array_almost_equal(mstats.normaltest(x), stats.normaltest(x))
assert_array_almost_equal(mstats.skewtest(x), stats.skewtest(x)) assert_array_almost_equal(mstats.skewtest(x), stats.skewtest(x))
assert_array_almost_equal(mstats.kurtosistest(x), assert_array_almost_equal(mstats.kurtosistest(x),
@ -534,7 +554,7 @@ class TestNormalitytests():
def test_axis_None(self): def test_axis_None(self):
# Test axis=None (equal to axis=0 for 1-D input) # Test axis=None (equal to axis=0 for 1-D input)
x = np.array((-2,-1,0,1,2,3)*4)**2 x = np.array((-2, -1, 0, 1, 2, 3) * 4) ** 2
assert_allclose(mstats.normaltest(x, axis=None), mstats.normaltest(x)) assert_allclose(mstats.normaltest(x, axis=None), mstats.normaltest(x))
assert_allclose(mstats.skewtest(x, axis=None), mstats.skewtest(x)) assert_allclose(mstats.skewtest(x, axis=None), mstats.skewtest(x))
assert_allclose(mstats.kurtosistest(x, axis=None), assert_allclose(mstats.kurtosistest(x, axis=None),
@ -542,7 +562,7 @@ class TestNormalitytests():
def test_maskedarray_input(self): def test_maskedarray_input(self):
# Add some masked values, test result doesn't change # Add some masked values, test result doesn't change
x = np.array((-2,-1,0,1,2,3)*4)**2 x = np.array((-2, -1, 0, 1, 2, 3) * 4) ** 2
xm = np.ma.array(np.r_[np.inf, x, 10], xm = np.ma.array(np.r_[np.inf, x, 10],
mask=np.r_[True, [False] * x.size, True]) mask=np.r_[True, [False] * x.size, True])
assert_allclose(mstats.normaltest(xm), stats.normaltest(x)) assert_allclose(mstats.normaltest(xm), stats.normaltest(x))
@ -550,7 +570,7 @@ class TestNormalitytests():
assert_allclose(mstats.kurtosistest(xm), stats.kurtosistest(x)) assert_allclose(mstats.kurtosistest(xm), stats.kurtosistest(x))
def test_nd_input(self): def test_nd_input(self):
x = np.array((-2,-1,0,1,2,3)*4)**2 x = np.array((-2, -1, 0, 1, 2, 3) * 4) ** 2
x_2d = np.vstack([x] * 2).T x_2d = np.vstack([x] * 2).T
for func in [mstats.normaltest, mstats.skewtest, mstats.kurtosistest]: for func in [mstats.normaltest, mstats.skewtest, mstats.kurtosistest]:
res_1d = func(x) res_1d = func(x)
@ -559,7 +579,7 @@ class TestNormalitytests():
assert_allclose(res_2d[1], [res_1d[1]] * 2) assert_allclose(res_2d[1], [res_1d[1]] * 2)
#TODO: for all ttest functions, add tests with masked array inputs # TODO: for all ttest functions, add tests with masked array inputs
class TestTtest_rel(): class TestTtest_rel():
def test_vs_nonmasked(self): def test_vs_nonmasked(self):

Loading…
Cancel
Save