Algorithms for calculating variance: Difference between revisions

:<math>\sigma^2 = \overline{(x^2)} - \bar x^2 = \frac{\sum_{i=1}^N x_i^2}{N} - \left(\frac{\sum_{i=1}^N x_i}{N}\right)^2</math>

<div style="margin-left: 35px; width: 600px">

{{framebox|blue}}

Because {{math|SumSq}} and {{math|(Sum×Sum)/''n''}} can be very similar numbers, [[Catastrophic cancellation|cancellation]] can lead to the [[precision (arithmetic)|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,<ref name="Einarsson2005">{{cite book|first=Bo|last=Einarsson|title=Accuracy and Reliability in Scientific Computing|url=https://books.google.com/books?id=8hrDV5EbrEsC|year=2005|publisher=SIAM|isbn=978-0-89871-584-2|page=47}}</ref><ref name="Chan1983">{{cite journal|url=http://cpsc.yale.edu/sites/default/files/files/tr222.pdf |archive-url=https://ghostarchive.org/archive/20221009/http://cpsc.yale.edu/sites/default/files/files/tr222.pdf |archive-date=2022-10-09 |url-status=live|first1=Tony F.|last1=Chan|author1-link=Tony F. Chan|first2=Gene H.|last2=Golub|author2-link=Gene H. Golub|first3=Randall J.|last3=LeVeque|title=Algorithms for computing the sample variance: Analysis and recommendations|journal=The American Statistician|volume=37|issue=3|pages=242–247|year=1983|jstor=2683386|doi=10.1080/00031305.1983.10483115}}</ref> and several alternate, numerically stable, algorithms have been proposed.<ref name=":1">{{Cite book|last1=Schubert|first1=Erich|last2=Gertz|first2=Michael|date=2018-07-09|title=Numerically stable parallel computation of (co-)variance|url=http://dl.acm.org/citation.cfm?id=3221269.3223036|publisher=ACM|pages=10|doi=10.1145/3221269.3223036|isbn=9781450365055|s2cid=49665540}}</ref> This is particularly bad if the standard deviation is small relative to the mean.

===Computing shifted data===

<syntaxhighlight lang="python">

 

 

 

 

 

 

==Two-pass algorithm==

This algorithm is numerically stable if ''n'' is small.<ref name="Einarsson2005"/><ref>{{cite book |last=Higham |first=Nicholas J. |url=https://epubs.siam.org/doi/book/10.1137/1.9780898718027 |title=Accuracy and Stability of Numerical Algorithms |publisher=Society for Industrial and Applied Mathematics |year=2002 |edition=2nd |publication-place=Philadelphia, PA |chapter=Problem 1.10 |doi= 10.1137/1.9780898718027|isbn=978-0-898715-21-7 |id=e{{ISBN|978-0-89871-802-7}}, 2002075848 }} Metadata also listed at [https://dl.acm.org/doi/10.5555/579525 ACM Digital Library].</ref> 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.

==Welford's online algorithm==

The following formulas can be used to update the [[mean]] and (estimated) variance of the sequence, for an additional element ''x''_''n''. Here, <math display="inline">\overline{x}_n = \frac{1}{n} \sum_{i=1}^n x_i </math> denotes the sample mean of the first ''n'' samples <math>(x_1,\dots,x_n)</math>, <math display="inline">\sigma^2_n = \frac{1}{n} \sum_{i=1}^n \left(x_i - \overline{x}_n \right)^2</math> their [[Variance#Biased_sample_variance|biased sample variance]], and <math display="inline">s^2_n = \frac{1}{n - 1} \sum_{i=1}^n \left(x_i - \overline{x}_n \right)^2</math> their [[Variance#Unbiased_sample_variance|unbiased sample variance]].

These formulas suffer from numerical instability {{cn|date=March 2023}}, as they repeatedly subtract a small number from a big number which scales with ''n''. A better quantity for updating is the sum of squares of differences from the current mean, <math display="inline">\sum_{i=1}^n (x_i - \bar x_n)^2</math>, here denoted <math>M_{2,n}</math>:

This algorithm was found by Welford,<ref>{{cite journal |first=B. P. |last=Welford |year=1962 |title=Note on a method for calculating corrected sums of squares and products |journal=[[Technometrics]] |volume=4 |issue=3 |pages=419–420 |jstor=1266577 |doi=10.2307/1266577}}</ref><ref>[[Donald E. Knuth]] (1998). ''[[The Art of Computer Programming]]'', volume 2: ''Seminumerical Algorithms'', 3rd edn., p.&nbsp;232. Boston: Addison-Wesley.</ref> and it has been thoroughly analyzed.<ref name="Chan1983" /><ref>{{cite journal |last=Ling |first=Robert F. |year=1974 |title=Comparison of Several Algorithms for Computing Sample Means and Variances |journal=Journal of the American Statistical Association |volume=69 |issue=348 |pages=859–866 |doi=10.2307/2286154|jstor=2286154 }}</ref> It is also common to denote <math>M_k = \bar x_k</math> and <math>S_k = M_{2,k}</math>.<ref>{{cite web |last=Cook |first=John D. |date=30 September 2022 |orig-date=1 November 2014 |title=Accurately computing sample variance |url=http://www.johndcook.com/standard_deviation.html |website=John D. Cook Consulting: Expert consulting in applied mathematics & data privacy}}</ref>

An example Python implementation for Welford's algorithm is given below.

def update(existing_aggregate, new_value):

    (count, mean, M2) = existing_aggregate

    count += 1

    delta = new_value - mean

    mean += delta / count

    delta2 = new_value - mean

     M2 += delta * delta2

    return (count, mean, M2)

def finalize(existing_aggregate):

    (count, mean, M2) = existing_aggregate

     if count < 2:

        return float("nan")

    else:

         (mean, variance, sample_variance) = (mean, M2 / count, M2 / (count - 1))

        return (mean, variance, sample_variance)

==Weighted incremental algorithm==

Chan's method for estimating the mean is numerically unstable when <math>n_A \approx n_B</math> and both are large, because the numerical error in <math>\delta = \bar x_B - \bar x_A</math> is not scaled down in the way that it is in the <math>n_B = 1</math> case. In such cases, prefer <math display="inline">\bar x_{AB} = \frac{n_A \bar x_A + n_B \bar x_B}{n_{AB}}</math>.

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)

This can be generalized to allow parallelization with [[Advanced Vector Extensions|AVX]], with [[Graphics processing unit|GPUs]], and [[Computer cluster|computer clusters]], and to covariance.<ref name=":1" />

Assume that all floating point operations use standard [[IEEE 754#Double-precision 64 bit|IEEE 754 double-precision]]{{Broken anchor|date=2025-06-11|bot=User:Cewbot/log/20201008/configuration|target_link=IEEE 754#Double-precision 64 bit|reason= }} 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.

Terriberry<ref>{{Cite web |last=Terriberry |first=Timothy B. |date=15 October 2008 |orig-date=9 December 2007 |title=Computing Higher-Order Moments Online |url=http://people.xiph.org/~tterribe/notes/homs.html |url-status=dead |archive-url=https://web.archive.org/web/20140423031833/http://people.xiph.org/~tterribe/notes/homs.html |archive-date=23 April 2014 |access-date=5 May 2008}}</ref> extends Chan's formulae to calculating the third and fourth [[central moment]]s, needed for example when estimating [[skewness]] and [[kurtosis]]:

:<math>

\begin{align}

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

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:

Very similar algorithms can be used to compute the [[covariance]].

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:

     for ix, iy in zip(data_x, data_y):

         Ex += ix - kx

         Ey += iy - ky

         Exy += (ix - kx) * (iy - ky)

     for i1, i2 in zip(data1, data2):

         a = i1 - mean1

     # Bessel's correction for sample variance

     sample_covar = C / (n - 1)