# Algorithms for calculating variance

Algorithms for calculating variance play a major role in computational statistics. A key difficulty 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

A formula for calculating the variance of an entire population of size N is:

${\displaystyle \sigma ^{2}={\bar {(x^{2})}}-{\bar {x}}^{2}=\displaystyle {\frac {\sum _{i=1}^{N}x_{i}^{2}-(\sum _{i=1}^{N}x_{i})^{2}/N}{N}}.\!}$

Using Bessel's correction to calculate an unbiased estimate of the population variance from a finite sample of n observations, the formula is:

${\displaystyle 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:

• Let n ← 0, Sum ← 0, SumSq ← 0
• For each datum x:
• nn + 1
• Sum ← Sum + x
• SumSq ← SumSq + x × x
• Var = (SumSq − (Sum × Sum) / n) / (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 SumSq and (Sum×Sum)/n can be very similar numbers, cancellation can lead to the precision of the result to be much less than the inherent precision of the floating-point arithmetic used to perform the computation. Thus this algorithm should not be used in practice.[1][2] This is particularly bad if the standard deviation is small relative to the mean. However, the algorithm can be improved by adopting the method of the assumed mean.

### Computing shifted data

We can use a property of the variance to avoid the catastrophic cancellation in this formula, namely the variance is invariant with respect to changes in a location parameter

${\displaystyle \operatorname {Var} (X-k)=\operatorname {Var} (X).}$

with ${\displaystyle k}$ any constant, which leads to the new formula

${\displaystyle s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-K)^{2}-(\sum _{i=1}^{n}(x_{i}-K))^{2}/n}{n-1}}.\!}$

the closer ${\displaystyle K}$ is to the mean value the more accurate the result will be, but just choosing a value inside the samples range will guarantee the desired stability. If the values ${\displaystyle (x_{i}-K)}$ are small then there are no problems with the sum of its squares, on the contrary, if they are large it necessarily means that the variance is large as well. In any case the second term in the formula is always smaller than the first one therefore no cancellation may occur.[2]

If we take just the first sample as ${\displaystyle K}$ the algorithm can be written in Python programming language as

def shifted_data_variance(data):
if len(data) == 0:
return 0
K = data[0]
n = 0
sum_ = 0
sum_sqr = 0
for x in data:
n = n + 1
sum_ += x - K
sum_sqr += (x - K) * (x - K)
variance = (sum_sqr - (sum_ * sum_)/n)/(n - 1)
# use n instead of (n-1) if want to compute the exact variance of the given data
# use (n-1) if data are samples of a larger population
return variance


this formula facilitates as well the incremental computation, that can be expressed as

K = 0
n = 0
ex = 0
ex2 = 0

if (n == 0):
K = x
n = n + 1
ex += x - K
ex2 += (x - K) * (x - K)

def remove_variable(x):
n = n - 1
ex -= (x - K)
ex2 -= (x - K) * (x - K)

def get_meanvalue():
return K + ex / n

def get_variance():
return (ex2 - (ex*ex)/n) / (n-1)


## Two-pass algorithm

An alternative approach, using a different formula for the variance, first computes the sample mean,

${\displaystyle {\bar {x}}=\displaystyle {\frac {\sum _{j=1}^{n}x_{j}}{n}}}$,

and then computes the sum of the squares of the differences from the mean,

${\displaystyle \mathrm {variance} =s^{2}=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}{n-1}}\!}$,

where s is the standard deviation. This is given by the following pseudocode:

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

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

mean = sum1 / n

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

variance = sum2 / (n - 1)
return variance


This algorithm is numerically stable if n is small.[1][3] However, the results of both of these simple algorithms ("Naïve" and "Two-pass") 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.

## Online algorithm

It is often useful to be able to compute the variance in a single pass, inspecting each value ${\displaystyle x_{i}}$ 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 xn. Here, xn denotes the sample mean of the first n samples (x1, ..., xn), s2n their sample variance, and σ2n their population variance.

${\displaystyle {\bar {x}}_{n}={\frac {(n-1)\,{\bar {x}}_{n-1}+x_{n}}{n}}={\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\!}$
${\displaystyle s_{n}^{2}={\frac {(n-2)}{(n-1)}}\,s_{n-1}^{2}+{\frac {(x_{n}-{\bar {x}}_{n-1})^{2}}{n}},\quad n>1}$
${\displaystyle \sigma _{n}^{2}={\frac {(n-1)\,\sigma _{n-1}^{2}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}{n}}.}$

These formulas suffer from numerical instability. A better quantity for updating is the sum of squares of differences from the current mean, ${\displaystyle \textstyle \sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})^{2}}$, here denoted ${\displaystyle M_{2,n}}$:

${\displaystyle M_{2,n}\!=M_{2,n-1}+(x_{n}-{\bar {x}}_{n-1})(x_{n}-{\bar {x}}_{n})}$
${\displaystyle s_{n}^{2}={\frac {M_{2,n}}{n-1}}}$
${\displaystyle \sigma _{n}^{2}={\frac {M_{2,n}}{n}}}$

A numerically stable algorithm for the sample variance is given below. It also computes the mean. This algorithm was found by Welford,[4][5] and it has been thoroughly analyzed.[6][7] It is also common to denote ${\displaystyle M_{k}={\bar {x}}_{k}}$ and ${\displaystyle S_{k}=M_{2,k}}$.[8]

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

for x in data:
n += 1
delta = x - mean
mean += delta/n
M2 += delta*(x - mean)

if n < 2:
return float('nan')
else:
return M2 / (n - 1)


This algorithm is much less prone to loss of precision due to catastrophic 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, one can first compute and subtract an estimate of the mean, and then use this algorithm on the residuals.

The parallel algorithm below illustrates how to merge multiple sets of statistics calculated online.

## Weighted incremental algorithm

The algorithm can be extended to handle unequal sample weights, replacing the simple counter n with the sum of weights seen so far. West (1979)[9] suggests this incremental algorithm:

def weighted_incremental_variance(dataWeightPairs):
sumweight = 0
mean = 0
M2 = 0

for x, weight in dataWeightPairs:  # Alternatively "for x, weight in zip(data, weights):"
temp = weight + sumweight
delta = x - mean
R = delta * weight / temp
mean += R
M2 += sumweight * delta * R  # Alternatively, "M2 = M2 + weight * delta * (x−mean)"
sumweight = temp

variance_n = M2/sumweight
variance = variance_n * len(dataWeightPairs)/(len(dataWeightPairs) - 1)


## Parallel algorithm

Chan et al.[10] note that the above "On-line" algorithm is a special case of an algorithm that works for any partition of the sample ${\displaystyle X}$ into sets ${\displaystyle X_{A}}$, ${\displaystyle X_{B}}$:

${\displaystyle \delta \!={\bar {x}}_{B}-{\bar {x}}_{A}}$
${\displaystyle {\bar {x}}_{X}={\bar {x}}_{A}+\delta \cdot {\frac {n_{B}}{n_{X}}}}$
${\displaystyle 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.

Chan's method for estimating the mean is numerically unstable when ${\displaystyle n_{A}\approx n_{B}}$ and both are large, because the numerical error in ${\displaystyle {\bar {x}}_{B}-{\bar {x}}_{A}}$ is not scaled down in the way that it is in the ${\displaystyle n_{B}=1}$ case. In such cases, prefer ${\displaystyle {\bar {x}}_{X}={\frac {n_{A}{\bar {x}}_{A}+n_{B}{\bar {x}}_{B}}{n_{A}+n_{B}}}}$.

def parallel_variance(avg_a, count_a, var_a, avg_b, count_b, var_b):
delta = avg_b - avg_a
m_a = var_a * (count_a - 1)
m_b = var_b * (count_b - 1)
M2 = m_a + m_b + delta ** 2 * count_a * count_b / (count_a + count_b)
return M2 / (count_a + count_b - 1)


## 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 "Naïve" algorithm and "Two-pass" algorithm 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. "Two-pass" algorithm computes this variance estimate correctly, but "Naïve" algorithm returns 29.333333333333332 instead of 30. While this loss of precision may be tolerable and viewed as a minor flaw of "Naïve" algorithm, 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 "Two-pass"" algorithm, but "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[11] extends Chan's formulae to calculating the third and fourth central moments, needed for example when estimating skewness and kurtosis:

${\displaystyle 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_{A}M_{2,B}-n_{B}M_{2,A}}{n_{X}}}}$
{\displaystyle {\begin{aligned}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_{A}M_{3,B}-n_{B}M_{3,A}}{n_{X}}}\\\end{aligned}}}

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

skewness: ${\displaystyle g_{1}={\frac {{\sqrt {n}}M_{3}}{M_{2}^{3/2}}},}$
kurtosis: ${\displaystyle g_{2}={\frac {nM_{4}}{M_{2}^{2}}}-3.}$

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

${\displaystyle \delta \!=x-m}$
${\displaystyle m'=m+{\frac {\delta }{n}}}$
${\displaystyle M_{2}'=M_{2}+\delta ^{2}{\frac {n-1}{n}}}$
${\displaystyle M_{3}'=M_{3}+\delta ^{3}{\frac {(n-1)(n-2)}{n^{2}}}-{\frac {3\delta M_{2}}{n}}}$
${\displaystyle 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 ${\displaystyle \delta /n}$, 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 = 0
mean = 0
M2 = 0
M3 = 0
M4 = 0

for x in data:
n1 = n
n = n + 1
delta = x - mean
delta_n = delta / n
delta_n2 = delta_n * delta_n
term1 = delta * delta_n * n1
mean = mean + delta_n
M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
M2 = M2 + term1

kurtosis = (n*M4) / (M2*M2) - 3
return kurtosis


Pébay[12] 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.

Choi and Sweetman [13] offer two alternative methods to compute the skewness and kurtosis, each of which can save substantial computer memory requirements and CPU time in certain applications. The first approach is to compute the statistical moments by separating the data into bins and then computing the moments from the geometry of the resulting histogram, which effectively becomes a one-pass algorithm for higher moments. One benefit is that the statistical moment calculations can be carried out to arbitrary accuracy such that the computations can be tuned to the precision of, e.g., the data storage format or the original measurement hardware. A relative histogram of a random variable can be constructed in the conventional way: the range of potential values is divided into bins and the number of occurrences within each bin are counted and plotted such that the area of each rectangle equals the portion of the sample values within that bin:

${\displaystyle H(x_{k})={\frac {h(x_{k})}{A}}}$

where ${\displaystyle h(x_{k})}$ and ${\displaystyle H(x_{k})}$ represent the frequency and the relative frequency at bin ${\displaystyle x_{k}}$ and ${\displaystyle A=\sum _{k=1}^{K}h(x_{k})\,\Delta x_{k}}$ is the total area of the histogram. After this normalization, the ${\displaystyle n}$ raw moments and central moments of ${\displaystyle x(t)}$ can be calculated from the relative histogram:

${\displaystyle m_{n}^{(h)}=\sum _{k=1}^{K}x_{k}^{n}\,H(x_{k})\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}x_{k}^{n}\,h(x_{k})\Delta x_{k}}$
${\displaystyle \theta _{n}^{(h)}=\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,H(x_{k})\Delta x_{k}={\frac {1}{A}}\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,h(x_{k})\Delta x_{k}}$

where the superscript ${\displaystyle ^{(h)}}$ indicates the moments are calculated from the histogram. For constant bin width ${\displaystyle \Delta x_{k}=\Delta x}$ these two expressions can be simplified using ${\displaystyle I=A/\Delta x}$:

${\displaystyle m_{n}^{(h)}={\frac {1}{I}}{\sum _{k=1}^{K}x_{k}^{n}\,h(x_{k})}}$
${\displaystyle \theta _{n}^{(h)}={\frac {1}{I}}{\sum _{k=1}^{K}{\Big (}x_{k}-m_{1}^{(h)}{\Big )}^{n}\,h(x_{k})}}$

The second approach from Choi and Sweetman [13] is an analytical methodology to combine statistical moments from individual segments of a time-history such that the resulting overall moments are those of the complete time-history. This methodology could be used for parallel computation of statistical moments with subsequent combination of those moments, or for combination of statistical moments computed at sequential times.

If ${\displaystyle Q}$ sets of statistical moments are known: ${\displaystyle (\gamma _{0,q},\mu _{q},\sigma _{q}^{2},\alpha _{3,q},\alpha _{4,q})\quad }$ for ${\displaystyle q=1,2,...,Q}$, then each ${\displaystyle \gamma _{n}}$ can be expressed in terms of the equivalent ${\displaystyle n}$ raw moments:

${\displaystyle \gamma _{n,q}=m_{n,q}\gamma _{0,q}\qquad \quad {\textrm {for}}\quad n=1,2,3,4\quad {\text{ and }}\quad q=1,2,\dots ,Q}$

where ${\displaystyle \gamma _{0,q}}$ is generally taken to be the duration of the ${\displaystyle q^{th}}$ time-history, or the number of points if ${\displaystyle \Delta t}$ is constant.

The benefit of expressing the statistical moments in terms of ${\displaystyle \gamma }$ is that the ${\displaystyle Q}$ sets can be combined by addition, and there is no upper limit on the value of ${\displaystyle Q}$.

${\displaystyle \gamma _{n,c}=\sum _{q=1}^{Q}\gamma _{n,q}\quad \quad {\textrm {for}}\quad n=0,1,2,3,4}$

where the subscript ${\displaystyle _{c}}$ represents the concatenated time-history or combined ${\displaystyle \gamma }$. These combined values of ${\displaystyle \gamma }$ can then be inversely transformed into raw moments representing the complete concatenated time-history

${\displaystyle m_{n,c}={\frac {\gamma _{n,c}}{\gamma _{0,c}}}\quad {\textrm {for}}\quad n=1,2,3,4}$

Known relationships between the raw moments (${\displaystyle m_{n}}$) and the central moments (${\displaystyle \theta _{n}=E[(x-\mu )^{n}])}$) are then used to compute the central moments of the concatenated time-history. Finally, the statistical moments of the concatenated history are computed from the central moments:

${\displaystyle \mu _{c}=m_{1,c}\ \ \ \ \ \sigma _{c}^{2}=\theta _{2,c}\ \ \ \ \ \alpha _{3,c}={\frac {\theta _{3,c}}{\sigma _{c}^{3}}}\ \ \ \ \ \alpha _{4,c}={\frac {\theta _{4,c}}{\sigma _{c}^{4}}}-3}$

## Covariance

Very similar algorithms can be used to compute the covariance. The naive algorithm is:

${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}x_{i}y_{i}-(\sum _{i=1}^{n}x_{i})(\sum _{i=1}^{n}y_{i})/n}{n}}.\!}$

For the algorithm above, one could use the following Python code:

def naive_covariance(data1, data2):
n = len(data1)
sum12 = 0
sum1 = sum(data1)
sum2 = sum(data2)

for i in range(n):
sum12 += data1[i]*data2[i]

covariance = (sum12 - sum1*sum2 / n) / n
return covariance


As for the variance, the covariance of two random variables is also shift-invariant, so given that ${\displaystyle K_{x}}$ and ${\displaystyle K_{y}}$ are whatever two constant values it can be written:

${\displaystyle \operatorname {Cov} (X,Y)=\operatorname {Cov} (X-k_{x},Y-k_{y})=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-K_{x})(y_{i}-K_{y})-(\sum _{i=1}^{n}(x_{i}-K_{x}))(\sum _{i=1}^{n}(y_{i}-K_{y}))/n}{n}}.\!}$

and again choosing a value inside the range of values will stabilize the formula against catastrophic cancellation as well as make it more robust against big sums. Taking the first value of each data set, the algorithm can be written as:

def shifted_data_covariance(dataX, dataY):
n = len(dataX)
if (n < 2):
return 0
Kx = dataX[0]
Ky = dataY[0]
Ex = 0
Ey = 0
Exy = 0
for i in range(n):
Ex += dataX[i] - Kx
Ey += dataY[i] - Ky
Exy += (dataX[i] - Kx) * (dataY[i] - Ky)
return (Exy - Ex * Ey / n) / n


The two-pass algorithm first computes the sample means, and then the covariance:

${\displaystyle {\bar {x}}=\displaystyle \sum _{i=1}^{n}x_{i}/n}$
${\displaystyle {\bar {y}}=\displaystyle \sum _{i=1}^{n}y_{i}/n}$
${\displaystyle \operatorname {Cov} (X,Y)=\displaystyle {\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{n}}.\!}$

The two-pass algorithm may be written as:

def two_pass_covariance(data1, data2):
n = len(data1)

mean1 = sum(data1) / n
mean2 = sum(data2) / n

covariance = 0

for i in range(n):
a = data1[i] - mean1
b = data2[i] - mean2
covariance += a*b / n

return covariance


A slightly more accurate compensated version performs the full naive algorithm on the residuals. The final sums ${\displaystyle \textstyle \sum x_{i}}$ and ${\displaystyle \textstyle \sum y_{i}}$ should be zero, but the second pass compensates for any small error.

A slight modification of the online algorithm for computing the variance yields an online algorithm for the covariance:

def online_covariance(data1, data2):
mean1 = mean2 = 0
M12 = 0
for x, y in zip(data1, data2):
n += 1
delta1 = (x - mean1) / n
mean1 += delta1
delta2 = (y - mean2) / n
mean2 += delta2
M12 += (n - 1) * delta1 * delta2 - M12 / n
return n / (n - 1) * M12


A stable one-pass algorithm exists, similar to the one above, that computes co-moment ${\displaystyle \textstyle C_{n}=\sum _{i=1}^{n}(x_{i}-{\bar {x}}_{n})(y_{i}-{\bar {y}}_{n})}$:

${\displaystyle {\bar {x}}_{n}={\bar {x}}_{n-1}+{\frac {x_{n}-{\bar {x}}_{n-1}}{n}}\!}$
${\displaystyle {\bar {y}}_{n}={\bar {y}}_{n-1}+{\frac {y_{n}-{\bar {y}}_{n-1}}{n}}\!}$
${\displaystyle C_{n}=C_{n-1}+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})=C_{n-1}+(y_{n}-{\bar {y}}_{n})(x_{n}-{\bar {x}}_{n-1})}$

The apparent asymmetry in that last equation is due to the fact that ${\displaystyle \textstyle (x_{n}-{\bar {x}}_{n})={\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})}$, so both update terms are equal to ${\displaystyle \textstyle {\frac {n-1}{n}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}$. Even greater accuracy can be achieved by first computing the means, then using the stable one-pass algorithm on the residuals.

Thus we can compute the covariance as

{\displaystyle {\begin{aligned}\operatorname {Cov} _{N}(X,Y)={\frac {C_{N}}{N}}&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(x_{n}-{\bar {x}}_{n})(y_{n}-{\bar {y}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+(y_{n}-{\bar {y}}_{n})(x_{n}-{\bar {x}}_{n-1})}{N}}\\&={\frac {\operatorname {Cov} _{N-1}(X,Y)\cdot (N-1)+{\frac {N-1}{N}}(x_{n}-{\bar {x}}_{n-1})(y_{n}-{\bar {y}}_{n-1})}{N}}.\end{aligned}}}

Likewise, there is a formula for combining the covariances of two sets that can be used to parallelize the computation:

${\displaystyle C_{X}=C_{A}+C_{B}+({\bar {x}}_{A}-{\bar {x}}_{B})({\bar {y}}_{A}-{\bar {y}}_{B})\cdot {\frac {n_{A}n_{B}}{n_{X}}}}$.

## References

1. ^ a b Bo Einarsson (1 August 2005). Accuracy and Reliability in Scientific Computing. SIAM. p. 47. ISBN 978-0-89871-584-2. Retrieved 17 February 2013.
2. ^ a b T.F.Chan, G.H. Golub and R.J. LeVeque (1983). ""Algorithms for computing the sample variance: Analysis and recommendations", The American Statistician, 37" (PDF): 242–247.
3. ^ Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10). SIAM.
4. ^ B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420.
5. ^ Donald E. Knuth (1998). The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
6. ^ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). Algorithms for Computing the Sample Variance: Analysis and Recommendations. The American Statistician 37, 242-247. http://www.jstor.org/stable/2683386
7. ^ Ling, Robert F. (1974). Comparison of Several Algorithms for Computing Sample Means and Variances. Journal of the American Statistical Association, Vol. 69, No. 348, 859-866. doi:10.2307/2286154
8. ^ http://www.johndcook.com/standard_deviation.html
9. ^ D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method
10. ^ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances." (PDF), Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University.
11. ^ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online
12. ^ Pébay, Philippe (2008), "Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments" (PDF), Technical Report SAND2008-6212, Sandia National Laboratories
13. ^ a b Choi, Muenkeun; Sweetman, Bert (2010), Efficient Calculation of Statistical Moments for Structural Health Monitoring (PDF)