# [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

trunk/scipy/stats/models/tests/test_scale.py
Modified:
trunk/scipy/stats/models/robust/scale.py
Log:

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 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
+
+
+    """
+    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:
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

===================================================================
--- 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):
+
+        X = W((40,10))
+        self.assertEquals(m.shape, (10,))
+
+        X = W((40,10,30))
+        self.assertEquals(m.shape, (10,30))
+
+        self.assertEquals(m.shape, (40,30))
+
+        self.assertEquals(m.shape, (40,10))
+
+        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()

```