I submitted this to the Sourceforge bug tracker, but wanted also to let this list know, as this is a potentially nasty bug. MLab.std() gives completely incorrect answers for multidimensional arrays when axis != 0.

...

...

foo array([[[ 1., 1., 1.], [ 2., 2., 2.], [ 3., 3., 3.]], [[ 1., 4., 4.], [ 2., 5., 5.], [ 3., 6., 6.]]]) std(foo) array([[ 0. , 2.12132034, 2.12132034], [ 0. , 2.12132034, 2.12132034], [ 0. , 2.12132034, 2.12132034]]) std(foo, 1) array([[ 0., 0., 0.], [ 0., 0., 0.]])

The following should fix the problem (but I haven't tested it extensively): def std(m,axis=0): """std(m,axis=0) returns the standard deviation along the given dimension of m. The result is unbiased with division by N-1. If m is of integer type returns a floating point answer. """ x = asarray(m) n = float(x.shape[axis]) x2 = mean(x * x, axis) x = mean(x, axis) return sqrt((x2 - x * x) * n /(n-1.0)) Jonathan Gilligan