the mean is numerically unstable when and both are large, because the numerical error in is not scaled down in the way that it is in the case. In such cases, prefer . def parallel_variance(n_a, avg_a, M2_a, n_b, avg_b, M2_b): n = n_a + n_b delta = avg_b - avg_a M2 = M2_a + M2_b + delta**2 * n_a * n_b / n var_ab = M2 / (n - 1) return var_ab This can be generalized to allow parallelization with AVX, with GPUs, and computer clusters, and to covariance. Example Assume that all floating point operations use standard IEEE 754 double-precision arithmetic. Consider the sample (4, 7, 13, 16) from an infinite population. Based on this sample, the estimated population mean is 10, and the unbiased estimate of population variance is 30. Both the naïve algorithm and two-pass algorithm compute these values correctly. Next consider the sample (, , , ), which gives rise to the same estimated variance as the first sample. The two-pass algorithm computes this variance estimate correctly, but the naïve algorithm returns 29.333333333333332 instead of 30. While this loss of precision may be tolerable and viewed as a minor flaw of the naïve algorithm, further increasing the offset makes the error catastrophic. Consider the sample (, , , ). Again the estimated population variance of 30 is computed correctly by the two-pass algorithm, but the naïve algorithm now computes it as −170.66666666666666. This is a serious problem with naïve algorithm and is due to catastrophic cancellation in the subtraction of two similar numbers at the final stage of the algorithm. Higher-order statistics Terriberry extends Chan's formulae to calculating the third and fourth central moments, needed for example when estimating skewness and kurtosis: Here the are again the sums of powers of differences from the mean , giving For the incremental case (i.e., ), this simplifies to: By preserving the value , only one division operation is needed and the higher-order statistics can thus be calculated for little incremental cost. An example of the online algorithm for kurtosis implemented as described is: def online_kurtosis(data): n = mean = M2 = M3 = M4 = 0 for x in data: n1 = n n = n + 1 delta = x - mean delta_n = delta / n delta_n2 = delta_n**2 term1 = delta * delta_n * n1 mean = mean + delta_n M4 = M4 + term1