[Numpy-discussion] Bug in MLab.std()

[Jonathan M. G. <jon...@va...>]

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