Algorithms_for_calculating

'''[[Algorithm]]s for calculating [[variance]]''' play a minor role in [[statistics|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.

== Algorithm I ==
The [[formula]] for calculating the variance of an entire [[statistical population|population]] of size ''n'' is:

:<math>\sigma^2 = \frac {\sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2/n}{n}. \!</math>

The formula for calculating an [[estimator bias|unbiased]] estimate of the population variance from a finite [[statistical sample|sample]] of ''n'' observations is:

:<math>s^2 = \frac {\sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2/n}{n-1}. \!</math>

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

n = 0
sum = 0
sum_sqr = 0

foreach x in data:
n = n + 1
sum = sum + x
sum_sqr = sum_sqr + x*x
end for

mean = sum/n
variance = (sum_sqr - sum*mean)/(n - 1)

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 <code>sum_sqr</code> and <code>sum * mean</code> can be very similar numbers, the [[precision (arithmetic)|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 variance is small relative to the sum of the numbers.

== Algorithm II ==
An alternate approach, using a different formula for the variance, is given by the following pseudocode:

n = 0
sum1 = 0
foreach x in data:
n = n + 1
sum1 = sum1 + x
end for
mean = sum1/n

sum2 = 0
foreach x in data:
sum2 = sum2 + (x - mean)^2
end for
variance = sum2/(n - 1)

This algorithm is often more numerically reliable than Algorithm I 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.

== Algorithm II (compensated) ==
The compensated-summation version of the algorithm above reads:

n = 0
sum1 = 0
foreach x in data:
n = n + 1
sum1 = sum1 + x
end for
mean = sum1/n

sum2 = 0
sumc = 0
foreach x in data:
sum2 = sum2 + (x - mean)^2
sumc = sumc + (x - mean)
end for
variance = (sum2 - sumc^2/n)/(n - 1)

== Algorithm III ==
The following formulas can be used to update the [[mean]] and (estimated) variance of the sequence, for an additional element <math>x_{\mathrm{new}}</math>. Here, ''m'' denotes the estimate of the population mean, ''s''2n-1 the estimate of the population variance, ''s''2n the estimate of the sample variance, and ''n'' the number of elements in the sequence before the addition.

:<math>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} \!</math>

:<math>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} \!</math>

:<math>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}.</math>

A numerically stable algorithm is given below. It also computes the mean.
This algorithm is due to Knuth,<ref>[[Donald E. Knuth]] (1998). ''[[The Art of Computer Programming]]'', volume 2: ''Seminumerical Algorithms'', 3rd edn., p. 232. Boston: Addison-Wesley.</ref>
who cites Welford.<ref>B. P. Welford (1962).[http://www.jstor.org/view/00401706/ap040015/04a00090/0 "Note on a method for calculating corrected sums of squares and products"]. ''[[Technometrics]]'' 4(3):419–420.</ref>

n = 0
mean = 0
S = 0

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

variance = S/(n - 1)

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.

== Algorithm IV ==
When the observations are weighted, West (1979) <ref>D. H. D. West (1979). ''[[Communications of the ACM]]'', 22, 9, 532-535: ''Updating Mean and Variance Estimates: An Improved Method''</ref> suggests this incremental algorithm:

n = 0
foreach x in the data:
if n=0 then
n = 1
mean = x
S = 0
sumweight = weight
else
n = n + 1
temp = weight + sumweight
S = S + sumweight*weight*(x-mean)^2 / temp
mean = mean + (x-mean)*weight / temp
sumweight = temp
end if
end for
Variance = S * sumweight * n / (n-1) // if sample is the population, omit n/(n-1)

== Example ==
Assume that all floating point operations use the standard [[IEEE 754#Double-precision 64 bit|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 (<math>10^8+4</math>, <math>10^8+7</math>, <math>10^8+13</math>, <math>10^8+16</math>), 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 (<math>10^9+4</math>, <math>10^9+7</math>, <math>10^9+13</math>, <math>10^9+16</math>). 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 because by definition the variance can never be negative.

== References ==

<references/>

== External links ==
* {{MathWorld|title=Sample Variance Computation|urlname=SampleVarianceComputation}}

[[Category:Statistical algorithms]]
[[Category:Articles with example pseudocode]]