FWIW, statistics.variance appears to do all its arithmetic using exact arbitrary-precision fractions, which in addition to being extremely slow will probably blow up your memory on reasonable-sized datasets.

If you want to check for 0 variance, you should use p2p

I see, the statistics approach is indeed too high a price to pay for numpy. I figured there might be a reformulated numerical algorithm where loss of precision does not cause the same problem but I might be wrong there.

I ran into it in a vectorised situation, where I would have noticed a bug in my code months ago if my .std(x) had been zero instead of small (due to things being infinite). I'm prepending a .ptp() check now to act correctly where ptp() is zero.

Perhaps I'm asking too much here. FWIW, I confirmed that Matlab also gives a nonzero answer (although not the same as numpy!).

I'm not aware of any general-purpose numerical algorithms that can give that kind of precision guarantee. I think the only ways you could guarantee exact 0 from a variance calculation are by using an infinite precision calculation, or by doing a preprocessing pass to explicitly check for all the values being equal. I guess the latter isn't even necessarily all that expensive in common cases where you can bail out after checking a few values, but it's extra complexity (esp. in the vectorized case) for a pretty specific and unusual case. I dunno.

I'm not aware of any general-purpose numerical algorithms that can give that kind of precision guarantee.

FWIW, it is not hard to show that Welford's algorithm [1] will give a variance that is exactly 0 when the input is constant.

Numpy uses a two-pass algorithm, and the underlying problem appears to be that numpy's calculation of the mean of x is off by one ULP:

In [258]: x = np.array([6715266981.538051]*10)

In [259]: np.mean(x)
Out[259]: 6715266981.5380497

In [260]: np.mean(x) == x[0]
Out[260]: False

[1] https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm

Ah, you're quite right -- I actually checked that one before posting, but misread the details :-). Do we know how it compares to our current algorithm in terms of speed or accuracy?

I don't know about the speed. John D. Cook blogged about the accuracy: https://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/

Related paper with accuracy and runtime comparisons of various methods: https://dl.acm.org/citation.cfm?id=3223036

@sk1p - thanks for the link! I think it would be great to use that one-pass algorithm. Perhaps ideally for a new gufunc that would calculate (weighted) mean and variance in one go.

Following up on this: a quick gufunc implementation by @0x0L of the algorithm referenced by @sk1p above correctly gives a variance of zero for the example on top.