Skip to content
This repository
Browse code

Merge branch 'pull-391-mstats' into master. Closes #1798.

Reviewed at #391
  • Loading branch information...
commit b69fe1860220e83bbdd52f7a0c4606c600c4ce4b 2 parents 785491c + 54f74c1
Ralf Gommers authored December 27, 2012
20  scipy/stats/mstats_basic.py
@@ -1434,7 +1434,12 @@ def skew(a, axis=0, bias=True):
1434 1434
     n = a.count(axis)
1435 1435
     m2 = moment(a, 2, axis)
1436 1436
     m3 = moment(a, 3, axis)
1437  
-    vals = ma.where(m2 == 0, 0, m3 / m2**1.5)
  1437
+    olderr = np.seterr(all='ignore')
  1438
+    try:
  1439
+        vals = ma.where(m2 == 0, 0, m3 / m2**1.5)
  1440
+    finally:
  1441
+        np.seterr(**olderr)
  1442
+
1438 1443
     if not bias:
1439 1444
         can_correct = (n > 2) & (m2 > 0)
1440 1445
         if can_correct.any():
@@ -1448,13 +1453,19 @@ def skew(a, axis=0, bias=True):
1448 1453
 
1449 1454
 def kurtosis(a, axis=0, fisher=True, bias=True):
1450 1455
     a, axis = _chk_asarray(a, axis)
1451  
-    n = a.count(axis)
1452 1456
     m2 = moment(a,2,axis)
1453 1457
     m4 = moment(a,4,axis)
1454  
-    vals = ma.where(m2 == 0, 0, m4/ m2**2.0)
  1458
+    olderr = np.seterr(all='ignore')
  1459
+    try:
  1460
+        vals = ma.where(m2 == 0, 0, m4 / m2**2.0)
  1461
+    finally:
  1462
+        np.seterr(**olderr)
  1463
+
1455 1464
     if not bias:
1456  
-        can_correct = (n > 3) & (m2 > 0)
  1465
+        n = a.count(axis)
  1466
+        can_correct = (n > 3) & (m2 is not ma.masked and m2 > 0)
1457 1467
         if can_correct.any():
  1468
+            n = np.extract(can_correct, n)
1458 1469
             m2 = np.extract(can_correct, m2)
1459 1470
             m4 = np.extract(can_correct, m4)
1460 1471
             nval = 1.0/(n-2)/(n-3)*((n*n-1.0)*m4/m2**2.0-3*(n-1)**2.0)
@@ -1465,6 +1476,7 @@ def kurtosis(a, axis=0, fisher=True, bias=True):
1465 1476
         return vals
1466 1477
 kurtosis.__doc__ = stats.kurtosis.__doc__
1467 1478
 
  1479
+
1468 1480
 def describe(a, axis=0):
1469 1481
     """
1470 1482
     Computes several descriptive statistics of the passed array.
10  scipy/stats/stats.py
@@ -1060,7 +1060,6 @@ def kurtosis(a, axis=0, fisher=True, bias=True):
1060 1060
         The kurtosis of values along an axis. If all values are equal,
1061 1061
         return -3 for Fisher's definition and 0 for Pearson's definition.
1062 1062
 
1063  
-
1064 1063
     References
1065 1064
     ----------
1066 1065
     [CRCProbStat2000]_ Section  2.2.25
@@ -1075,7 +1074,12 @@ def kurtosis(a, axis=0, fisher=True, bias=True):
1075 1074
     m2 = moment(a,2,axis)
1076 1075
     m4 = moment(a,4,axis)
1077 1076
     zero = (m2 == 0)
1078  
-    vals = np.where(zero, 0, m4/ m2**2.0)
  1077
+    olderr = np.seterr(all='ignore')
  1078
+    try:
  1079
+        vals = np.where(zero, 0, m4 / m2**2.0)
  1080
+    finally:
  1081
+        np.seterr(**olderr)
  1082
+
1079 1083
     if not bias:
1080 1084
         can_correct = (n > 3) & (m2 > 0)
1081 1085
         if can_correct.any():
@@ -3034,7 +3038,7 @@ def ttest_ind(a, b, axis=0, equal_var=True):
3034 3038
         population variance [2]_.
3035 3039
 
3036 3040
         .. versionadded:: 0.11.0
3037  
-  
  3041
+
3038 3042
     Returns
3039 3043
     -------
3040 3044
     t : float or array
83  scipy/stats/tests/test_mstats_basic.py
@@ -9,9 +9,10 @@
9 9
 from numpy.ma import masked, nomask
10 10
 
11 11
 import scipy.stats.mstats as mstats
  12
+from scipy import stats
12 13
 from numpy.testing import TestCase, run_module_suite
13 14
 from numpy.ma.testutils import assert_equal, assert_almost_equal, \
14  
-    assert_array_almost_equal, assert_
  15
+    assert_array_almost_equal, assert_array_almost_equal_nulp, assert_
15 16
 
16 17
 
17 18
 class TestMquantiles(TestCase):
@@ -36,7 +37,6 @@ def test_mquantiles_limit_keyword(self):
36 37
         assert_almost_equal(quants, desired)
37 38
 
38 39
 
39  
-
40 40
 class TestGMean(TestCase):
41 41
     def test_1D(self):
42 42
         a = (1,2,3,4)
@@ -72,6 +72,7 @@ def test_2D(self):
72 72
                             np.power(1*4,1./2.)))
73 73
         assert_array_almost_equal(actual, desired, decimal=14)
74 74
 
  75
+
75 76
 class TestHMean(TestCase):
76 77
     def test_1D(self):
77 78
         a = (1,2,3,4)
@@ -277,21 +278,30 @@ def test_winsorization(self):
277 278
 
278 279
 
279 280
 class TestMoments(TestCase):
280  
-    """
281  
-        Comparison numbers are found using R v.1.5.1
282  
-        note that length(testcase) = 4
283  
-        testmathworks comes from documentation for the
284  
-        Statistics Toolbox for Matlab and can be found at both
285  
-        http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml
286  
-        http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml
287  
-        Note that both test cases came from here.
288  
-    """
  281
+    # Comparison numbers are found using R v.1.5.1
  282
+    # note that length(testcase) = 4
  283
+    # testmathworks comes from documentation for the
  284
+    # Statistics Toolbox for Matlab and can be found at both
  285
+    # http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml
  286
+    # http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml
  287
+    # Note that both test cases came from here.
289 288
     testcase = [1,2,3,4]
290 289
     testmathworks = ma.fix_invalid([1.165 , 0.6268, 0.0751, 0.3516, -0.6965,
291 290
                                     np.nan])
  291
+    testcase_2d = ma.array(
  292
+    np.array([[ 0.05245846,  0.50344235,  0.86589117,  0.36936353,  0.46961149],
  293
+           [ 0.11574073,  0.31299969,  0.45925772,  0.72618805,  0.75194407],
  294
+           [ 0.67696689,  0.91878127,  0.09769044,  0.04645137,  0.37615733],
  295
+           [ 0.05903624,  0.29908861,  0.34088298,  0.66216337,  0.83160998],
  296
+           [ 0.64619526,  0.94894632,  0.27855892,  0.0706151 ,  0.39962917]]),
  297
+    mask = np.array([[ True, False, False,  True, False],
  298
+           [ True,  True,  True, False,  True],
  299
+           [False, False, False, False, False],
  300
+           [True, True, True, True, True],
  301
+           [False, False,  True, False, False]], dtype=np.bool))
  302
+
292 303
     def test_moment(self):
293  
-        """
294  
-        mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))"""
  304
+        # mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))
295 305
         y = mstats.moment(self.testcase,1)
296 306
         assert_almost_equal(y,0.0,10)
297 307
         y = mstats.moment(self.testcase,2)
@@ -300,17 +310,16 @@ def test_moment(self):
300 310
         assert_almost_equal(y,0.0)
301 311
         y = mstats.moment(self.testcase,4)
302 312
         assert_almost_equal(y,2.5625)
  313
+
303 314
     def test_variation(self):
304  
-        """variation = samplestd/mean """
  315
+        #variation = samplestd/mean
305 316
 ##        y = stats.variation(self.shoes[0])
306 317
 ##        assert_almost_equal(y,21.8770668)
307 318
         y = mstats.variation(self.testcase)
308 319
         assert_almost_equal(y,0.44721359549996, 10)
309 320
 
310 321
     def test_skewness(self):
311  
-        """
312  
-            sum((testmathworks-mean(testmathworks,axis=0))**3,axis=0)/((sqrt(var(testmathworks)*4/5))**3)/5
313  
-        """
  322
+        # sum((testmathworks-mean(testmathworks,axis=0))**3,axis=0)/((sqrt(var(testmathworks)*4/5))**3)/5
314 323
         y = mstats.skew(self.testmathworks)
315 324
         assert_almost_equal(y,-0.29322304336607,10)
316 325
         y = mstats.skew(self.testmathworks,bias=0)
@@ -319,26 +328,46 @@ def test_skewness(self):
319 328
         assert_almost_equal(y,0.0,10)
320 329
 
321 330
     def test_kurtosis(self):
322  
-        """
323  
-            sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
324  
-            sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
325  
-            Set flags for axis = 0 and
326  
-            fisher=0 (Pearson's definition of kurtosis for compatibility with Matlab)
327  
-        """
  331
+        #    sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4
  332
+        #    sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5
  333
+        #    Set flags for axis = 0 and
  334
+        #    fisher=0 (Pearson's definition of kurtosis for compatibility with Matlab)
328 335
         y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1)
329 336
         assert_almost_equal(y, 2.1658856802973,10)
330 337
         # Note that MATLAB has confusing docs for the following case
331 338
         #  kurtosis(x,0) gives an unbiased estimate of Pearson's skewness
332 339
         #  kurtosis(x)  gives a biased estimate of Fisher's skewness (Pearson-3)
333 340
         #  The MATLAB docs imply that both should give Fisher's
334  
-        y = mstats.kurtosis(self.testmathworks,fisher=0,bias=0)
  341
+        y = mstats.kurtosis(self.testmathworks,fisher=0, bias=0)
335 342
         assert_almost_equal(y, 3.663542721189047,10)
336 343
         y = mstats.kurtosis(self.testcase,0,0)
337 344
         assert_almost_equal(y,1.64)
338  
-    #
  345
+
  346
+        # test that kurtosis works on multidimensional masked arrays
  347
+        correct_2d = ma.array(np.array([-1.5, -3., -1.47247052385,  0.,
  348
+                                        -1.26979517952]),
  349
+                              mask=np.array([False, False, False,  True,
  350
+                                             False], dtype=np.bool))
  351
+        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1),
  352
+                                  correct_2d)
  353
+        for i, row in enumerate(self.testcase_2d):
  354
+            assert_almost_equal(mstats.kurtosis(row), correct_2d[i])
  355
+
  356
+        correct_2d_bias_corrected = ma.array(
  357
+            np.array([-1.5, -3., -1.88988209538,  0., -0.5234638463918877]),
  358
+            mask=np.array([False, False, False,  True, False], dtype=np.bool))
  359
+        assert_array_almost_equal(mstats.kurtosis(self.testcase_2d, 1,
  360
+                                                  bias=False),
  361
+                                  correct_2d_bias_corrected)
  362
+        for i, row in enumerate(self.testcase_2d):
  363
+            assert_almost_equal(mstats.kurtosis(row, bias=False),
  364
+                                correct_2d_bias_corrected[i])
  365
+
  366
+        # Check consistency between stats and mstats implementations
  367
+        assert_array_almost_equal_nulp(mstats.kurtosis(self.testcase_2d[2, :]),
  368
+                                       stats.kurtosis(self.testcase_2d[2, :]))
  369
+
339 370
     def test_mode(self):
340  
-        "Tests the mode"
341  
-        #
342 371
         a1 = [0,0,0,1,1,1,2,3,3,3,3,4,5,6,7]
343 372
         a2 = np.reshape(a1, (3,5))
344 373
         ma1 = ma.masked_where(ma.array(a1)>2,a1)

0 notes on commit b69fe18

Please sign in to comment.
Something went wrong with that request. Please try again.