The way I compute the variance above is, to my knowledge, the standard algorithm implemented by most software packages (apparently including MATLAB). (The single pass computation you mention is subject to catastrophic cancellation, and thus pretty terrible unless the mean of your data is very small relative to the variance.) However, the "real" implementation in Julia standard library is indeed a bit better: it performs pairwise summation, which has O(log n) error growth instead of O(n) error growth at negligible performance cost (see http://en.wikipedia.org/wiki/Pairwise_summation).