Algorithms for calculating variance: Wikis

Note: Many of our articles have direct quotes from sources you can cite, within the Wikipedia article! This article doesn't yet, but we're working on it! See more info or our list of citable articles.

Encyclopedia

Algorithms for calculating variance play a major role in statistical computing. A key problem in the design of good algorithms for this problem is that formulas for the variance may involve sums of squares, which can lead to numerical instability as well as to arithmetic overflow when dealing with large values.

Naïve algorithm

The formula for calculating the variance of an entire population of size n is:

$\sigma^2 = \displaystyle\frac {\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2/n}{n}. \!$

The formula for calculating an unbiased estimate of the population variance from a finite sample of n observations is:

$s^2 = \displaystyle\frac {\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2/n}{n-1}. \!$

Therefore a naive algorithm to calculate the estimated variance is given by the following:

def naive_variance(data):
n = 0
Sum = 0
Sum_sqr = 0

for x in data:
n = n + 1
Sum = Sum + x
Sum_sqr = Sum_sqr + x*x

mean = Sum/n
variance = (Sum_sqr - Sum*mean)/(n - 1)
return variance


This algorithm can easily be adapted to compute the variance of a finite population: simply divide by n instead of n − 1 on the last line.

Because sum_sqr and sum * mean can be very similar numbers, the precision of the result can be much less than the inherent precision of the floating-point arithmetic used to perform the computation. This is particularly bad if the standard deviation is small relative to the mean.

Two-pass algorithm

An alternate approach, using a different formula for the variance,

$s^2 = \displaystyle\frac {\sum_{i=1}^n (x_i^2 - \sum_{j=1}^n x_j/n)^2}{n-1}. \!$

is given by the following pseudocode:

def two_pass_variance(data):
n = 0
sum1 = 0

for x in data:
n = n + 1
sum1 = sum1 + x

mean = sum1/n

sum2 = 0
for x in data:
sum2 = sum2 + (x - mean)**2

variance = sum2/(n - 1)
return variance


This algorithm is often more numerically reliable than the naïve algorithm for large sets of data, although it can be worse if much of the data is very close to but not precisely equal to the mean and some are quite far away from it.

The results of both of these simple algorithms (I and II) can depend inordinately on the ordering of the data and can give poor results for very large data sets due to repeated roundoff error in the accumulation of the sums. Techniques such as compensated summation can be used to combat this error to a degree.

Compensated variant

The compensated-summation version of the algorithm above reads:

def compensated_variance(data):
n = 0
sum1 = 0
for x in data:
n = n + 1
sum1 = sum1 + x
mean = sum1/n

sum2 = 0
sumc = 0
for x in data:
sum2 = sum2 + (x - mean)**2
sumc = sumc + (x - mean)
variance = (sum2 - sumc**2/n)/(n - 1)
return variance


On-line algorithm

It is often useful to be able to compute the variance in a single pass, inspecting each value xi only once; for example, when the data are being collected without enough storage to keep all the values, or when costs of memory access dominate those of computation. For such an online algorithm, a recurrence relation is required between quantities from which the required statistics can be calculated in a numerically stable fashion.

The following formulas can be used to update the mean and (estimated) variance of the sequence, for an additional element xnew. Here, m denotes the estimate of the population mean(using the sample mean), s2n-1 the estimate of the sample variance, s2n the estimate of the population variance, and n the number of elements in the sequence before the addition.

$m_{\mathrm{new}} = \frac{n \; m_{\mathrm{old}} + x_{\mathrm{new}}}{n+1} = m_{\mathrm{old}} + \frac{x_{\mathrm{new}} - m_{\mathrm{old}}}{n+1} \!$
$s^2_{\mathrm{n-1, new}} = \frac{(n-1) \; s^2_{\mathrm{n-1, old}} + (x_{\mathrm{new}} - m_{\mathrm{new}}) \, (x_{\mathrm{new}} - m_{\mathrm{old}})}{n} \; \; \, \, \; \, \mathrm{n>0} \!$
$s^2_{\mathrm{n, new}} = \frac{n \; s^2_{\mathrm{n, old}} + (x_{\mathrm{new}} - m_{\mathrm{new}}) \, (x_{\mathrm{new}} - m_{\mathrm{old}})}{n+1}.$

It turns out that a more suitable quantity for updating is the sum of squares of differences from the (current) mean, $\sum_{i=1}^n (x_i - m)^2$, here denoted M2:

$M_\mathrm{2,new}\! = M_\mathrm{2,old} + (x_\mathrm{new} - m_\mathrm{old})(x_\mathrm{new} - m_\mathrm{new})$
$s^2_\mathrm{n} = \frac{M_2}{n}$
$s^2_\mathrm{n-1} = \frac{M_2}{n-1}$

A numerically stable algorithm is given below. It also computes the mean. This algorithm is due to Knuth,[1] who cites Welford.[2]

def online_variance(data):
n = 0
mean = 0
M2 = 0

for x in data:
n = n + 1
delta = x - mean
mean = mean + delta/n
M2 = M2 + delta*(x - mean)  # This expression uses the new value of mean

variance_n = M2/n
variance = M2/(n - 1)
return variance


This algorithm is much less prone to loss of precision due to massive cancellation, but might not be as efficient because of the division operation inside the loop. For a particularly robust two-pass algorithm for computing the variance, first compute and subtract an estimate of the mean, and then use this algorithm on the residuals.

A slightly more convenient form allows one to calculate the standard deviation without having to explicitly calculate the new mean. If n is the number of elements in the sequence after the addition of the new element, then one has

$s^2_{\mathrm{n-1, new}} = \frac{(n-2)s^2_{\mathrm{n-1, old}}+\frac{n-1}{n}\left(x_\text{new}-m_\text{old}\right)^2}{n-1}$

Weighted incremental algorithm

When the observations are weighted, West (1979)[3] suggests this incremental algorithm:

def weighted_incremental_variance(dataWeightPairs):
n = 0
mean = 0
S = 0
sumweight = 0
for x, weight in dataWeightPairs: # Alternately "for x in zip(data, weight):"
n = n + 1
temp = weight + sumweight
Q = x - mean
R = Q * weight / temp
S = S + sumweight * Q * R
mean = mean + R
sumweight = temp
Variance = S * n / ((n-1) * sumweight)  # if sample is the population, omit n/(n-1)
return Variance


Parallel algorithm

Chan et al.[4] note that the above on-line algorithm III is a special case of an algorithm that works for any partition of the sample X into sets XA, XB:

$\delta\! = m^B - m^A$
$m^X = m^A + \delta\cdot\frac{N^B}{N^X}$
$M_2^X = M_2^A + M_2^B + \delta^2\cdot\frac{N^A N^B}{N^X}$.

This may be useful when, for example, multiple processing units may be assigned to discrete parts of the input.

Higher-order statistics

Terriberry[5] extends Chan's formulae to calculating the third and fourth central moments, needed for example when estimating skewness and kurtosis:

$M_3^X = M_3^A + M_3^B + \delta^3\frac{N^A N^B (N^A - N^B)}{(N^X)^2} + 3\delta\frac{N^AM_2^B - N^BM_2^A}{N^X}$
\begin{align} M_4^X = M_4^A + M_4^B & + \delta^4\frac{N^A N^B \left((N^A)^2 - N^A N^B + (N^B)^2\right)}{(N^X)^3} \ & + 6\delta^2\frac{(N^A)^2 M_2^B + (N^B)^2 M_2^A}{(N^X)^2} + 4\delta\frac{N^AM_3^B - N^BM_3^A}{N^X} \ \end{align}

Here the Mk are again the sums of powers of differences from the mean $\Sigma(x - \overline{x})^k$, giving

skewness: $g_1 = \frac{\sqrt{n} M_3}{M_2^{3/2}},$
kurtosis: $g_2 = \frac{n M_4}{M_2^2}.$

For the incremental case (i.e., B = {x}), this simplifies to:

$\delta\! = x - m$
$m' = m + \frac{\delta}{n}$
$M_2' = M_2 + \delta^2 \frac{ n-1}{n}$
$M_3' = M_3 + \delta^3 \frac{ (n - 1) (n - 2)}{n^2} - \frac{3\delta M_2}{n}$
$M_4' = M_4 + \frac{\delta^4 (n - 1) (n^2 - 3n + 3)}{n^3} + \frac{6\delta^2 M_2}{n^2} - \frac{4\delta M_3}{n}$

By preserving the value δ / n, only one division operation is needed and the higher-order statistics can thus be calculated for little incremental cost.

Pébay[6] further extends these results to arbitrary-order central moments, for the incremental and the pairwise cases. One can also find there similar formulas for covariance.

Example

Assume that all floating point operations use the 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 Algorithm I and Algorithm II compute these values correctly. Next consider the sample (108 + 4, 108 + 7, 108 + 13, 108 + 16), which gives rise to the same estimated variance as the first sample. Algorithm II computes this variance estimate correctly, but Algorithm I returns 29.333333333333332 instead of 30. While this loss of precision may be tolerable and viewed as a minor flaw of Algorithm I, it is easy to find data that reveal a major flaw in the naive algorithm: Take the sample to be (109 + 4, 109 + 7, 109 + 13, 109 + 16). Again the estimated population variance of 30 is computed correctly by Algorithm II, but the naive algorithm now computes it as −170.66666666666666. This is a serious problem with Algorithm I and is due to catastrophic cancellation in the subtraction of two similar numbers at the final stage of the algorithm.

References

1. ^ Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
2. ^ B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420.
3. ^ D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method
4. ^ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances.", Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University  .
5. ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online
6. ^ Pébay, Philippe (2008), "Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments", Technical Report SAND2008-6212, Sandia National Laboratories