[Scipy-svn] r3443 - in trunk/scipy/stats/models: robust tests

scipy-svn@scip... scipy-svn@scip...
Wed Oct 17 13:28:35 CDT 2007


Author: jonathan.taylor
Date: 2007-10-17 13:28:32 -0500 (Wed, 17 Oct 2007)
New Revision: 3443

Added:
   trunk/scipy/stats/models/tests/test_scale.py
Modified:
   trunk/scipy/stats/models/robust/scale.py
Log:
added axis option to MAD estimate, and Huber estimate of scale

added test for scipy.models.robust.scale



Modified: trunk/scipy/stats/models/robust/scale.py
===================================================================
--- trunk/scipy/stats/models/robust/scale.py	2007-10-17 08:13:12 UTC (rev 3442)
+++ trunk/scipy/stats/models/robust/scale.py	2007-10-17 18:28:32 UTC (rev 3443)
@@ -1,19 +1,41 @@
 import numpy as N
-from scipy import median
-from scipy.stats import norm
+from scipy.stats import norm, median
 
-def MAD(a, c=0.6745):
+def unsqueeze(data, axis, oldshape):
     """
-    Median Absolute Deviation along first axis of an array:
+    unsqueeze a collapsed array
 
+    >>> from numpy import mean
+    >>> from numpy.random import standard_normal
+    >>> x = standard_normal((3,4,5))
+    >>> m = mean(x, axis=1)
+    >>> m.shape
+    (3, 5)
+    >>> unsqueeze(m, 1, x.shape)
+    >>> m.shape
+    (3, 1, 5)
+    >>>                       
+    """
+    
+    newshape = list(oldshape)
+    newshape[axis] = 1
+    data.shape = newshape
+    
+
+def MAD(a, c=0.6745, axis=0):
+    """
+    Median Absolute Deviation along given axis of an array:
+
     median(abs(a - median(a))) / c
 
     """
 
     a = N.asarray(a, N.float64)
-    d = N.multiply.outer(median(a), N.ones(a.shape[1:]))
-    return median(N.fabs(a - d) / c)
+    d = median(a, axis=axis)
+    unsqueeze(d, axis, a.shape)
 
+    return median(N.fabs(a - d) / c, axis=axis)
+
 class Huber:
     """
     Huber's proposal 2 for estimating scale.
@@ -29,19 +51,20 @@
     gamma = tmp + c**2 * (1 - tmp) - 2 * c * norm.pdf(c)
     del(tmp)
     
-    niter = 10
+    niter = 30
 
-    def __call__(self, a, mu=None, scale=None):
+    def __call__(self, a, mu=None, scale=None, axis=0):
         """
         Compute Huber\'s proposal 2 estimate of scale, using an optional
         initial value of scale and an optional estimate of mu. If mu
         is supplied, it is not reestimated.
         """
 
+        self.axis = axis
         self.a = N.asarray(a, N.float64)
         if mu is None:
             self.n = self.a.shape[0] - 1
-            self.mu = N.multiply.outer(median(self.a), N.ones(self.a.shape[1:]))
+            self.mu = median(self.a, axis=axis)
             self.est_mu = True
         else:
             self.n = self.a.shape[0]
@@ -49,14 +72,18 @@
             self.est_mu = False
 
         if scale is None:
-            self.scale = MAD(self.a)**2
+            self.scale = MAD(self.a, axis=self.axis)**2
         else:
             self.scale = scale
 
+        unsqueeze(self.scale, self.axis, self.a.shape)
+        unsqueeze(self.mu, self.axis, self.a.shape)
+
         for donothing in self:
             pass
 
-        self.s = N.sqrt(self.scale)
+        self.s = N.squeeze(N.sqrt(self.scale))
+        del(self.scale); del(self.mu); del(self.a)
         return self.s
 
     def __iter__(self):
@@ -67,16 +94,16 @@
         a = self.a
         subset = self.subset(a)
         if self.est_mu:
-            mu = (subset * a).sum() / a.shape[0]
+            mu = N.sum(subset * a + (1 - Huber.c) * subset, axis=self.axis) / a.shape[self.axis]
         else:
             mu = self.mu
+        unsqueeze(mu, self.axis, self.a.shape)
             
-        scale = N.sum(subset * (a - mu)**2, axis=0) / (self.n * Huber.gamma - N.sum(1. - subset, axis=0) * Huber.c**2)
+        scale = N.sum(subset * (a - mu)**2, axis=self.axis) / (self.n * Huber.gamma - N.sum(1. - subset, axis=self.axis) * Huber.c**2)
 
         self.iter += 1
 
-        if (N.fabs(N.sqrt(scale) - N.sqrt(self.scale)) <= N.sqrt(self.scale) * Huber.tol and
-            N.fabs(mu - self.mu) <= N.sqrt(self.scale) * Huber.tol):
+        if N.alltrue(N.less_equal(N.fabs(N.sqrt(scale) - N.sqrt(self.scale)), N.sqrt(self.scale) * Huber.tol)) and N.alltrue(N.less_equal(N.fabs(mu - self.mu), N.sqrt(self.scale) * Huber.tol)):
             self.scale = scale
             self.mu = mu
             raise StopIteration
@@ -84,6 +111,8 @@
             self.scale = scale
             self.mu = mu
 
+        unsqueeze(self.scale, self.axis, self.a.shape)
+
         if self.iter >= self.niter:
             raise StopIteration
 

Added: trunk/scipy/stats/models/tests/test_scale.py
===================================================================
--- trunk/scipy/stats/models/tests/test_scale.py	2007-10-17 08:13:12 UTC (rev 3442)
+++ trunk/scipy/stats/models/tests/test_scale.py	2007-10-17 18:28:32 UTC (rev 3443)
@@ -0,0 +1,53 @@
+"""
+Test functions for models.robust.scale
+"""
+
+import numpy.random as R
+from numpy.testing import NumpyTest, NumpyTestCase
+
+import scipy.stats.models.robust.scale as scale
+
+W = R.standard_normal
+
+class TestScale(NumpyTestCase):
+
+    def test_MAD(self):
+        X = W((40,10))
+	m = scale.MAD(X)
+        self.assertEquals(m.shape, (10,))
+
+    def test_MADaxes(self):
+        X = W((40,10,30))
+	m = scale.MAD(X, axis=0)
+        self.assertEquals(m.shape, (10,30))
+
+	m = scale.MAD(X, axis=1)
+        self.assertEquals(m.shape, (40,30))
+
+	m = scale.MAD(X, axis=2)
+        self.assertEquals(m.shape, (40,10))
+
+	m = scale.MAD(X, axis=-1)
+        self.assertEquals(m.shape, (40,10))
+
+    def test_huber(self):
+        X = W((40,10))
+	m = scale.huber(X)
+        self.assertEquals(m.shape, (10,))
+
+    def test_huberaxes(self):
+        X = W((40,10,30))
+	m = scale.huber(X, axis=0)
+        self.assertEquals(m.shape, (10,30))
+
+	m = scale.huber(X, axis=1)
+        self.assertEquals(m.shape, (40,30))
+
+	m = scale.huber(X, axis=2)
+        self.assertEquals(m.shape, (40,10))
+
+	m = scale.huber(X, axis=-1)
+        self.assertEquals(m.shape, (40,10))
+
+if __name__ == "__main__":
+    NumpyTest().run()



More information about the Scipy-svn mailing list