[Scipy-svn] r3481 - in trunk/scipy/stats: . tests

[scipy-svn at scipy.org]

Author: cookedm
Date: 2007-11-01 09:25:15 -0500 (Thu, 01 Nov 2007)
New Revision: 3481

Modified:
   trunk/scipy/stats/stats.py
   trunk/scipy/stats/tests/test_stats.py
Log:
Fix #526 -- scipy.stats.gmean cannot handle large numbers
 - also minor fixes in stats/tests/test_stats.py


Modified: trunk/scipy/stats/stats.py
===================================================================
--- trunk/scipy/stats/stats.py	2007-11-01 11:56:59 UTC (rev 3480)
+++ trunk/scipy/stats/stats.py	2007-11-01 14:25:15 UTC (rev 3481)
@@ -348,8 +348,9 @@
 
     Parameters
     ----------
-    a : array
+    a : array of positive values
     axis : int or None
+    zero_sub : value to substitute for zero values. Default is 0.
 
     Returns
     -------
@@ -357,12 +358,11 @@
     all values in the array if axis==None.
     """
     a, axis = _chk_asarray(a, axis)
-    size = a.shape[axis]
-    prod = np.product(a, axis)
-    return np.power(prod, 1./size)
+    log_a = np.log(a)
+    return np.exp(log_a.mean(axis=axis))
 
 
-def hmean(a, axis=0):
+def hmean(a, axis=0, zero_sub=0):
     """Calculates the harmonic mean of the values in the passed array.
 
     That is:  n / (1/x1 + 1/x2 + ... + 1/xn)

Modified: trunk/scipy/stats/tests/test_stats.py
===================================================================
--- trunk/scipy/stats/tests/test_stats.py	2007-11-01 11:56:59 UTC (rev 3480)
+++ trunk/scipy/stats/tests/test_stats.py	2007-11-01 14:25:15 UTC (rev 3481)
@@ -500,18 +500,19 @@
         a = (1,2,3,4)
         actual= stats.gmean(a)
         desired = power(1*2*3*4,1./4.)
-        assert_almost_equal(desired,actual,decimal=14)
+        assert_almost_equal(actual, desired,decimal=14)
 
         desired1 = stats.gmean(a,axis=-1)
-        assert_almost_equal(desired1,actual,decimal=14)
+        assert_almost_equal(actual, desired1, decimal=14)
+
     def check_1D_array(self):
         a = array((1,2,3,4), float32)
         actual= stats.gmean(a)
         desired = power(1*2*3*4,1./4.)
-        assert_almost_equal(desired,actual,decimal=7)
+        assert_almost_equal(actual, desired, decimal=7)
 
         desired1 = stats.gmean(a,axis=-1)
-        assert_almost_equal(desired1,actual,decimal=7)
+        assert_almost_equal(actual, desired1, decimal=7)
 
     def check_2D_array_default(self):
         a = array(((1,2,3,4),
@@ -519,10 +520,10 @@
                    (1,2,3,4)))
         actual= stats.gmean(a)
         desired = array((1,2,3,4))
-        assert_array_almost_equal(desired,actual,decimal=14)
+        assert_array_almost_equal(actual, desired, decimal=14)
 
         desired1 = stats.gmean(a,axis=0)
-        assert_array_almost_equal(desired1,actual,decimal=14)
+        assert_array_almost_equal(actual, desired1, decimal=14)
 
     def check_2D_array_dim1(self):
         a = array(((1,2,3,4),
@@ -531,25 +532,30 @@
         actual= stats.gmean(a, axis=1)
         v = power(1*2*3*4,1./4.)
         desired = array((v,v,v))
-        assert_array_almost_equal(desired,actual,decimal=14)
+        assert_array_almost_equal(actual, desired, decimal=14)
 
+    def check_large_values(self):
+        a = array([1e100, 1e200, 1e300])
+        actual = stats.gmean(a)
+        assert_approx_equal(actual, 1e200, significant=14)
+
 class TestHMean(NumpyTestCase):
     def check_1D_list(self):
         a = (1,2,3,4)
         actual= stats.hmean(a)
         desired =  4. / (1./1 + 1./2 + 1./3 + 1./4)
-        assert_almost_equal(desired,actual,decimal=14)
+        assert_almost_equal(actual, desired, decimal=14)
 
         desired1 = stats.hmean(array(a),axis=-1)
-        assert_almost_equal(desired1,actual,decimal=14)
+        assert_almost_equal(actual, desired1, decimal=14)
     def check_1D_array(self):
         a = array((1,2,3,4), float64)
         actual= stats.hmean(a)
         desired =  4. / (1./1 + 1./2 + 1./3 + 1./4)
-        assert_almost_equal(desired,actual,decimal=14)
+        assert_almost_equal(actual, desired, decimal=14)
 
         desired1 = stats.hmean(a,axis=-1)
-        assert_almost_equal(desired1,actual,decimal=14)
+        assert_almost_equal(actual, desired1, decimal=14)
 
     def check_2D_array_default(self):
         a = array(((1,2,3,4),
@@ -557,10 +563,10 @@
                    (1,2,3,4)))
         actual = stats.hmean(a)
         desired = array((1.,2.,3.,4.))
-        assert_array_almost_equal(desired,actual,decimal=14)
+        assert_array_almost_equal(actual, desired, decimal=14)
 
         actual1 = stats.hmean(a,axis=0)
-        assert_array_almost_equal(desired,actual1,decimal=14)
+        assert_array_almost_equal(actual1, desired, decimal=14)
 
     def check_2D_array_dim1(self):
         a = array(((1,2,3,4),
@@ -570,7 +576,7 @@
         v = 4. / (1./1 + 1./2 + 1./3 + 1./4)
         desired1 = array((v,v,v))
         actual1 = stats.hmean(a, axis=1)
-        assert_array_almost_equal(desired1,actual1,decimal=14)
+        assert_array_almost_equal(actual1, desired1, decimal=14)
 
 
 class TestMean(NumpyTestCase):
@@ -592,16 +598,18 @@
         a = [[1.0, 2.0, 3.0],
              [2.0, 4.0, 6.0],
              [8.0, 12.0, 7.0]]
-        A = array(a,'d')
-        N1,N2 = (3,3)
-        mn1 = zeros(N2,'d')
+        A = array(a)
+        N1, N2 = (3, 3)
+        mn1 = zeros(N2, dtype=float)
         for k in range(N1):
             mn1 += A[k,:] / N1
-        allclose(stats.mean(a),mn1,rtol=1e-13,atol=1e-13)
-        mn2 = zeros(N1,'d')
+        assert_almost_equal(stats.mean(a, axis=0), mn1, decimal=13)
+        assert_almost_equal(stats.mean(a), mn1, decimal=13)
+        mn2 = zeros(N1, dtype=float)
         for k in range(N2):
-            mn2 += A[:,k] / N2
-        allclose(stats.mean(a,axis=0),mn2,rtol=1e-13,atol=1e-13)
+            mn2 += A[:,k]
+        mn2 /= N2
+        assert_almost_equal(stats.mean(a, axis=1), mn2, decimal=13)
 
     def check_ravel(self):
         a = rand(5,3,5)