Algorithms for calculating variance

Related on Lably

Through Algorithm, Variance and Sampling (statistics)

Lab · Class 8–12Algorithms & Sorting Lab · Class 11Algorithm Visualizer Lab · Class 11Decision Trees & Random Forests Lab · Class 10–11Reinforcement Learning Basics

Naïve algorithm

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

σ^{2} = \bar{(x - \bar{x})^{2}} = \bar{(x^{2})} - {\bar{x}}^{2} = \frac{\sum_{i = 1}^{N} x_{i}^{2}}{N} - {(\frac{\sum_{i = 1}^{N} x_{i}}{N})}^{2}

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

s^{2} = (\frac{\sum_{i = 1}^{n} x_{i}^{2}}{n} - {(\frac{\sum_{i = 1}^{n} x_{i}}{n})}^{2}) \cdot \frac{n}{n - 1} .

Therefore, a naïve algorithm to calculate the estimated variance is given by the following:

Let n ← 0, Sum ← 0, SumSq ← 0
For each datum x:
- n ← n + 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, and several alternate, numerically stable, algorithms have been proposed. This is particularly bad if the standard deviation is small relative to the mean.

6 more chapters are free with a Lably account.

Computing shifted data

The variance is invariant with respect to changes in a location parameter, a property which can be used to avoid the catastrophic cancellation in this formula.

Var (X - K) = Var (X) .

with $K$ any constant, which leads to the new formula

σ^{2} = \frac{\sum_{i = 1}^{n} (x_{i} - K)^{2} - (\sum_{i = 1}^{n} (x_{i} - K))^{2} / n}{n - 1} .

The closer $K$ is to the mean value, the more accurate the result will be, but just choosing any value inside the sample range will provide the desired stability. If the values $(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.

The algorithm can be written in Python programming language as

<span id="mwTQ"></span><span class="k" id="mwTg">def</span><span class="w" id="mwTw"> </span><span class="nf" id="mwUA">shifted_data_variance</span><span class="p" id="mwUQ">(</span><span class="n" id="mwUg">data</span><span class="p" id="mwUw">,</span> <span class="n" id="mwVA">K</span><span class="p" id="mwVQ">):</span> <span class="k" id="mwVg">if</span> <span class="nb" id="mwVw">len</span><span class="p" id="mwWA">(</span><span class="n" id="mwWQ">data</span><span class="p" id="mwWg">)</span> <span class="o" id="mwWw"><</span> <span class="mi" id="mwXA">2</span><span class="p" id="mwXQ">:</span> <span class="k" id="mwXg">return</span> <span class="mf" id="mwXw">0.0</span> <span class="n" id="mwYA">n</span> <span class="o" id="mwYQ">=</span> <span class="n" id="mwYg">Ex</span> <span class="o" id="mwYw">=</span> <span class="n" id="mwZA">Ex2</span> <span class="o" id="mwZQ">=</span> <span class="mf" id="mwZg">0.0</span> <span class="k" id="mwZw">for</span> <span class="n" id="mwaA">x</span> <span class="ow" id="mwaQ">in</span> <span class="n" id="mwag">data</span><span class="p" id="mwaw">:</span> <span class="n" id="mwbA">n</span> <span class="o" id="mwbQ">+=</span> <span class="mi" id="mwbg">1</span> <span class="n" id="mwbw">Ex</span> <span class="o" id="mwcA">+=</span> <span class="n" id="mwcQ">x</span> <span class="o" id="mwcg">-</span> <span class="n" id="mwcw">K</span> <span class="n" id="mwdA">Ex2</span> <span class="o" id="mwdQ">+=</span> <span class="p" id="mwdg">(</span><span class="n" id="mwdw">x</span> <span class="o" id="mweA">-</span> <span class="n" id="mweQ">K</span><span class="p" id="mweg">)</span> <span class="o" id="mwew">**</span> <span class="mi" id="mwfA">2</span> <span class="n" id="mwfQ">variance</span> <span class="o" id="mwfg">=</span> <span class="p" id="mwfw">(</span><span class="n" id="mwgA">Ex2</span> <span class="o" id="mwgQ">-</span> <span class="n" id="mwgg">Ex</span><span class="o" id="mwgw">**</span><span class="mi" id="mwhA">2</span> <span class="o" id="mwhQ">/</span> <span class="n" id="mwhg">n</span><span class="p" id="mwhw">)</span> <span class="o" id="mwiA">/</span> <span class="p" id="mwiQ">(</span><span class="n" id="mwig">n</span> <span class="o" id="mwiw">-</span> <span class="mi" id="mwjA">1</span><span class="p" id="mwjQ">)</span> <span class="c1" id="mwjg"># use n instead of (n-1) if want to compute the exact variance of the given data</span> <span class="c1" id="mwjw"># use (n-1) if data are samples of a larger population</span> <span class="k" id="mwkA">return</span> <span class="n" id="mwkQ">variance</span>

Two-pass algorithm

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

\bar{x} = \frac{\sum_{j = 1}^{n} x_{j}}{n},

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

sample variance = s^{2} = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}}{n - 1},

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

<span id="mwnw"></span><span class="k" id="mwoA">def</span><span class="w" id="mwoQ"> </span><span class="nf" id="mwog">two_pass_variance</span><span class="p" id="mwow">(</span><span class="n" id="mwpA">data</span><span class="p" id="mwpQ">):</span> <span class="n" id="mwpg">n</span> <span class="o" id="mwpw">=</span> <span class="nb" id="mwqA">len</span><span class="p" id="mwqQ">(</span><span class="n" id="mwqg">data</span><span class="p" id="mwqw">)</span> <span class="n" id="mwrA">mean</span> <span class="o" id="mwrQ">=</span> <span class="nb" id="mwrg">sum</span><span class="p" id="mwrw">(</span><span class="n" id="mwsA">data</span><span class="p" id="mwsQ">)</span> <span class="o" id="mwsg">/</span> <span class="n" id="mwsw">n</span> <span class="n" id="mwtA">variance</span> <span class="o" id="mwtQ">=</span> <span class="nb" id="mwtg">sum</span><span class="p" id="mwtw">((</span><span class="n" id="mwuA">x</span> <span class="o" id="mwuQ">-</span> <span class="n" id="mwug">mean</span><span class="p" id="mwuw">)</span> <span class="o" id="mwvA">**</span> <span class="mi" id="mwvQ">2</span> <span class="k" id="mwvg">for</span> <span class="n" id="mwvw">x</span> <span class="ow" id="mwwA">in</span> <span class="n" id="mwwQ">data</span><span class="p" id="mwwg">)</span> <span class="o" id="mwww">/</span> <span class="p" id="mwxA">(</span><span class="n" id="mwxQ">n</span> <span class="o" id="mwxg">-</span> <span class="mi" id="mwxw">1</span><span class="p" id="mwyA">)</span> <span class="k" id="mwyQ">return</span> <span class="n" id="mwyg">variance</span>

This piece of code, if run on CPython 3.12 and newer, is always numerically stable. This is because these versions use Neumaier's compensated summation scheme for the sum() function, making it resistant to repeated roundoff error. Many numerically oriented languages provide similar constructs.

This algorithm, if implemented with a naive addition sum = 0; for x in data: sum += x; mean = sum / n, would only be numerically stable if n is small because of repeated round-off when adding a small number to a much, much larger one.

Incremental algorithms

It is often useful to be able to compute the variance in a single pass, inspecting each value $x_{i}$ only once; for example, when the data is 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.

Incremental shifted data

Our shifted-data algorithm from before includes the simple updating of data in a loop, which is naturally incremental. It can be expressed as:

<span id="mw7Q"></span><span class="k" id="mw7g">class</span><span class="w" id="mw7w"> </span><span class="nc" id="mw8A">ShiftDataVariance</span><span class="p" id="mw8Q">:</span> <span class="k" id="mw8g">def</span><span class="w" id="mw8w"> </span><span class="fm" id="mw9A">__init__</span><span class="p" id="mw9Q">(</span><span class="bp" id="mw9g">self</span><span class="p" id="mw9w">):</span> <span class="bp" id="mw-A">self</span><span class="o" id="mw-Q">.</span><span class="n" id="mw-g">K</span> <span class="o" id="mw-w">=</span> <span class="mf" id="mw_A">0.0</span> <span class="bp" id="mw_Q">self</span><span class="o" id="mw_g">.</span><span class="n" id="mw_w">n</span> <span class="o" id="mwAQA">=</span> <span class="mi" id="mwAQE">0</span> <span class="bp" id="mwAQI">self</span><span class="o" id="mwAQM">.</span><span class="n" id="mwAQQ">Ex</span> <span class="o" id="mwAQU">=</span> <span class="mf" id="mwAQY">0.0</span> <span class="bp" id="mwAQc">self</span><span class="o" id="mwAQg">.</span><span class="n" id="mwAQk">Ex2</span> <span class="o" id="mwAQo">=</span> <span class="mf" id="mwAQs">0.0</span> <span class="k" id="mwAQw">def</span><span class="w" id="mwAQ0"> </span><span class="nf" id="mwAQ4">add_variable</span><span class="p" id="mwAQ8">(</span><span class="bp" id="mwARA">self</span><span class="p" id="mwARE">,</span> <span class="n" id="mwARI">x</span><span class="p" id="mwARM">:</span> <span class="nb" id="mwARQ">float</span><span class="p" id="mwARU">):</span> <span class="k" id="mwARY">if</span> <span class="bp" id="mwARc">self</span><span class="o" id="mwARg">.</span><span class="n" id="mwARk">n</span> <span class="o" id="mwARo">==</span> <span class="mi" id="mwARs">0</span><span class="p" id="mwARw">:</span> <span class="bp" id="mwAR0">self</span><span class="o" id="mwAR4">.</span><span class="n" id="mwAR8">K</span> <span class="o" id="mwASA">=</span> <span class="n" id="mwASE">x</span> <span class="bp" id="mwASI">self</span><span class="o" id="mwASM">.</span><span class="n" id="mwASQ">n</span> <span class="o" id="mwASU">+=</span> <span class="mi" id="mwASY">1</span> <span class="bp" id="mwASc">self</span><span class="o" id="mwASg">.</span><span class="n" id="mwASk">Ex</span> <span class="o" id="mwASo">+=</span> <span class="n" id="mwASs">x</span> <span class="o" id="mwASw">-</span> <span class="bp" id="mwAS0">self</span><span class="o" id="mwAS4">.</span><span class="n" id="mwAS8">K</span> <span class="bp" id="mwATA">self</span><span class="o" id="mwATE">.</span><span class="n" id="mwATI">Ex2</span> <span class="o" id="mwATM">+=</span> <span class="p" id="mwATQ">(</span><span class="n" id="mwATU">x</span> <span class="o" id="mwATY">-</span> <span class="bp" id="mwATc">self</span><span class="o" id="mwATg">.</span><span class="n" id="mwATk">K</span><span class="p" id="mwATo">)</span> <span class="o" id="mwATs">**</span> <span class="mi" id="mwATw">2</span> <span class="k" id="mwAT0">def</span><span class="w" id="mwAT4"> </span><span class="nf" id="mwAT8">remove_variable</span><span class="p" id="mwAUA">(</span><span class="bp" id="mwAUE">self</span><span class="p" id="mwAUI">,</span> <span class="n" id="mwAUM">x</span><span class="p" id="mwAUQ">:</span> <span class="nb" id="mwAUU">float</span><span class="p" id="mwAUY">):</span> <span class="bp" id="mwAUc">self</span><span class="o" id="mwAUg">.</span><span class="n" id="mwAUk">n</span> <span class="o" id="mwAUo">-=</span> <span class="mi" id="mwAUs">1</span> <span class="bp" id="mwAUw">self</span><span class="o" id="mwAU0">.</span><span class="n" id="mwAU4">Ex</span> <span class="o" id="mwAU8">-=</span> <span class="n" id="mwAVA">x</span> <span class="o" id="mwAVE">-</span> <span class="bp" id="mwAVI">self</span><span class="o" id="mwAVM">.</span><span class="n" id="mwAVQ">K</span> <span class="bp" id="mwAVU">self</span><span class="o" id="mwAVY">.</span><span class="n" id="mwAVc">Ex2</span> <span class="o" id="mwAVg">-=</span> <span class="p" id="mwAVk">(</span><span class="n" id="mwAVo">x</span> <span class="o" id="mwAVs">-</span> <span class="bp" id="mwAVw">self</span><span class="o" id="mwAV0">.</span><span class="n" id="mwAV4">K</span><span class="p" id="mwAV8">)</span> <span class="o" id="mwAWA">**</span> <span class="mi" id="mwAWE">2</span> <span class="k" id="mwAWI">def</span><span class="w" id="mwAWM"> </span><span class="nf" id="mwAWQ">add_variables</span><span class="p" id="mwAWU">(</span><span class="bp" id="mwAWY">self</span><span class="p" id="mwAWc">,</span> <span class="n" id="mwAWg">xs</span><span class="p" id="mwAWk">:</span> <span class="nb" id="mwAWo">list</span><span class="p" id="mwAWs">[</span><span class="nb" id="mwAWw">float</span><span class="p" id="mwAW0">]):</span> <span class="c1" id="mwAW4"># Python uses Neumaier's algorithm for the builtin sum() function,</span> <span class="c1" id="mwAW8"># which is more accurate than a simple loop.</span> <span class="k" id="mwAXA">if</span> <span class="bp" id="mwAXE">self</span><span class="o" id="mwAXI">.</span><span class="n" id="mwAXM">n</span> <span class="o" id="mwAXQ">==</span> <span class="mi" id="mwAXU">0</span> <span class="ow" id="mwAXY">and</span> <span class="n" id="mwAXc">xs</span><span class="p" id="mwAXg">:</span> <span class="bp" id="mwAXk">self</span><span class="o" id="mwAXo">.</span><span class="n" id="mwAXs">K</span> <span class="o" id="mwAXw">=</span> <span class="n" id="mwAX0">xs</span><span class="p" id="mwAX4">[</span><span class="mi" id="mwAX8">0</span><span class="p" id="mwAYA">]</span> <span class="bp" id="mwAYE">self</span><span class="o" id="mwAYI">.</span><span class="n" id="mwAYM">n</span> <span class="o" id="mwAYQ">+=</span> <span class="nb" id="mwAYU">len</span><span class="p" id="mwAYY">(</span><span class="n" id="mwAYc">xs</span><span class="p" id="mwAYg">)</span> <span class="bp" id="mwAYk">self</span><span class="o" id="mwAYo">.</span><span class="n" id="mwAYs">Ex</span> <span class="o" id="mwAYw">+=</span> <span class="nb" id="mwAY0">sum</span><span class="p" id="mwAY4">(</span><span class="n" id="mwAY8">x</span> <span class="o" id="mwAZA">-</span> <span class="bp" id="mwAZE">self</span><span class="o" id="mwAZI">.</span><span class="n" id="mwAZM">K</span> <span class="k" id="mwAZQ">for</span> <span class="n" id="mwAZU">x</span> <span class="ow" id="mwAZY">in</span> <span class="n" id="mwAZc">xs</span><span class="p" id="mwAZg">)</span> <span class="bp" id="mwAZk">self</span><span class="o" id="mwAZo">.</span><span class="n" id="mwAZs">Ex2</span> <span class="o" id="mwAZw">+=</span> <span class="nb" id="mwAZ0">sum</span><span class="p" id="mwAZ4">((</span><span class="n" id="mwAZ8">x</span> <span class="o" id="mwAaA">-</span> <span class="bp" id="mwAaE">self</span><span class="o" id="mwAaI">.</span><span class="n" id="mwAaM">K</span><span class="p" id="mwAaQ">)</span> <span class="o" id="mwAaU">**</span> <span class="mi" id="mwAaY">2</span> <span class="k" id="mwAac">for</span> <span class="n" id="mwAag">x</span> <span class="ow" id="mwAak">in</span> <span class="n" id="mwAao">xs</span><span class="p" id="mwAas">)</span> <span class="k" id="mwAaw">def</span><span class="w" id="mwAa0"> </span><span class="nf" id="mwAa4">get_mean</span><span class="p" id="mwAa8">(</span><span class="bp" id="mwAbA">self</span><span class="p" id="mwAbE">)</span> <span class="o" id="mwAbI">-></span> <span class="nb" id="mwAbM">float</span><span class="p" id="mwAbQ">:</span> <span class="k" id="mwAbU">return</span> <span class="bp" id="mwAbY">self</span><span class="o" id="mwAbc">.</span><span class="n" id="mwAbg">K</span> <span class="o" id="mwAbk">+</span> <span class="bp" id="mwAbo">self</span><span class="o" id="mwAbs">.</span><span class="n" id="mwAbw">Ex</span> <span class="o" id="mwAb0">/</span> <span class="bp" id="mwAb4">self</span><span class="o" id="mwAb8">.</span><span class="n" id="mwAcA">n</span> <span class="k" id="mwAcE">def</span><span class="w" id="mwAcI"> </span><span class="nf" id="mwAcM">get_variance</span><span class="p" id="mwAcQ">(</span><span class="bp" id="mwAcU">self</span><span class="p" id="mwAcY">)</span> <span class="o" id="mwAcc">-></span> <span class="nb" id="mwAcg">float</span><span class="p" id="mwAck">:</span> <span class="k" id="mwAco">return</span> <span class="p" id="mwAcs">(</span><span class="bp" id="mwAcw">self</span><span class="o" id="mwAc0">.</span><span class="n" id="mwAc4">Ex2</span> <span class="o" id="mwAc8">-</span> <span class="bp" id="mwAdA">self</span><span class="o" id="mwAdE">.</span><span class="n" id="mwAdI">Ex</span><span class="o" id="mwAdM">**</span><span class="mi" id="mwAdQ">2</span> <span class="o" id="mwAdU">/</span> <span class="bp" id="mwAdY">self</span><span class="o" id="mwAdc">.</span><span class="n" id="mwAdg">n</span><span class="p" id="mwAdk">)</span> <span class="o" id="mwAdo">/</span> <span class="p" id="mwAds">(</span><span class="bp" id="mwAdw">self</span><span class="o" id="mwAd0">.</span><span class="n" id="mwAd4">n</span> <span class="o" id="mwAd8">-</span> <span class="mi" id="mwAeA">1</span><span class="p" id="mwAeE">)</span>

Welford's online algorithm

Analogous to the two-pass algorithm above, it remains desirable to use the actual mean of the data instead of just picking the first element. The following formulas can be used to update the mean and (estimated) variance of the sequence, for an additional element x_n. Here, ${\bar{x}}_{n} = \frac{1}{n} \sum_{i = 1}^{n} x_{i}$ denotes the sample mean of the first n samples $(x_{1}, \dots, x_{n})$ , $σ_{n}^{2} = \frac{1}{n} \sum_{i = 1}^{n} {(x_{i} - {\bar{x}}_{n})}^{2}$ their biased sample variance, and $s_{n}^{2} = \frac{1}{n - 1} \sum_{i = 1}^{n} {(x_{i} - {\bar{x}}_{n})}^{2}$ their unbiased sample variance.

{\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}

σ_{n}^{2} = \frac{(n - 1) σ_{n - 1}^{2} + (x_{n} - {\bar{x}}_{n - 1}) (x_{n} - {\bar{x}}_{n})}{n} = σ_{n - 1}^{2} + \frac{(x_{n} - {\bar{x}}_{n - 1}) (x_{n} - {\bar{x}}_{n}) - σ_{n - 1}^{2}}{n} .

s_{n}^{2} = \frac{n - 2}{n - 1} s_{n - 1}^{2} + \frac{(x_{n} - {\bar{x}}_{n - 1})^{2}}{n} = s_{n - 1}^{2} + \frac{(x_{n} - {\bar{x}}_{n - 1})^{2}}{n} - \frac{s_{n - 1}^{2}}{n - 1}, n > 1

These formulas suffer from numerical instability, 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, $\sum_{i = 1}^{n} (x_{i} - {\bar{x}}_{n})^{2}$ , here denoted $M_{2, n}$ :

\begin{matrix} M_{2, n} & = M_{2, n - 1} + (x_{n} - {\bar{x}}_{n - 1}) (x_{n} - {\bar{x}}_{n}) \\ σ_{n}^{2} & = \frac{M_{2, n}}{n} \\ s_{n}^{2} & = \frac{M_{2, n}}{n - 1} \end{matrix}

The following algorithm was found by Welford, and it has been thoroughly analyzed. It is also common to denote $M_{k} = {\bar{x}}_{k}$ and $S_{k} = M_{2, k}$ . An example Python implementation for Welford's algorithm is given below, using the same framework as the above "shifted data" algorithm:

<span id="mwAhk"></span><span class="k" id="mwAho">class</span><span class="w" id="mwAhs"> </span><span class="nc" id="mwAhw">WelfordVariance</span><span class="p" id="mwAh0">:</span> <span class="k" id="mwAh4">def</span><span class="w" id="mwAh8"> </span><span class="fm" id="mwAiA">__init__</span><span class="p" id="mwAiE">(</span><span class="bp" id="mwAiI">self</span><span class="p" id="mwAiM">):</span> <span class="c1" id="mwAiQ"># Comparison to ShiftDataVariance:</span> <span class="bp" id="mwAiU">self</span><span class="o" id="mwAiY">.</span><span class="n" id="mwAic">mean</span> <span class="o" id="mwAig">=</span> <span class="mf" id="mwAik">0.0</span> <span class="c1" id="mwAio"># = K + Ex / n</span> <span class="bp" id="mwAis">self</span><span class="o" id="mwAiw">.</span><span class="n" id="mwAi0">count</span> <span class="o" id="mwAi4">=</span> <span class="mi" id="mwAi8">0</span> <span class="c1" id="mwAjA"># = n</span> <span class="bp" id="mwAjE">self</span><span class="o" id="mwAjI">.</span><span class="n" id="mwAjM">M2</span> <span class="o" id="mwAjQ">=</span> <span class="mf" id="mwAjU">0.0</span> <span class="c1" id="mwAjY"># = Ex2 - (Ex)^2 / n</span> <span class="k" id="mwAjc">def</span><span class="w" id="mwAjg"> </span><span class="nf" id="mwAjk">add_variable</span><span class="p" id="mwAjo">(</span><span class="bp" id="mwAjs">self</span><span class="p" id="mwAjw">,</span> <span class="n" id="mwAj0">x</span><span class="p" id="mwAj4">:</span> <span class="nb" id="mwAj8">float</span><span class="p" id="mwAkA">):</span> <span class="bp" id="mwAkE">self</span><span class="o" id="mwAkI">.</span><span class="n" id="mwAkM">count</span> <span class="o" id="mwAkQ">+=</span> <span class="mi" id="mwAkU">1</span> <span class="n" id="mwAkY">old_mean</span> <span class="o" id="mwAkc">=</span> <span class="bp" id="mwAkg">self</span><span class="o" id="mwAkk">.</span><span class="n" id="mwAko">mean</span> <span class="bp" id="mwAks">self</span><span class="o" id="mwAkw">.</span><span class="n" id="mwAk0">mean</span> <span class="o" id="mwAk4">+=</span> <span class="p" id="mwAk8">(</span><span class="n" id="mwAlA">x</span> <span class="o" id="mwAlE">-</span> <span class="bp" id="mwAlI">self</span><span class="o" id="mwAlM">.</span><span class="n" id="mwAlQ">mean</span><span class="p" id="mwAlU">)</span> <span class="o" id="mwAlY">/</span> <span class="bp" id="mwAlc">self</span><span class="o" id="mwAlg">.</span><span class="n" id="mwAlk">count</span> <span class="bp" id="mwAlo">self</span><span class="o" id="mwAls">.</span><span class="n" id="mwAlw">M2</span> <span class="o" id="mwAl0">+=</span> <span class="p" id="mwAl4">(</span><span class="n" id="mwAl8">x</span> <span class="o" id="mwAmA">-</span> <span class="n" id="mwAmE">old_mean</span><span class="p" id="mwAmI">)</span> <span class="o" id="mwAmM">*</span> <span class="p" id="mwAmQ">(</span><span class="n" id="mwAmU">x</span> <span class="o" id="mwAmY">-</span> <span class="bp" id="mwAmc">self</span><span class="o" id="mwAmg">.</span><span class="n" id="mwAmk">mean</span><span class="p" id="mwAmo">)</span> <span class="k" id="mwAms">def</span><span class="w" id="mwAmw"> </span><span class="nf" id="mwAm0">remove_variable</span><span class="p" id="mwAm4">(</span><span class="bp" id="mwAm8">self</span><span class="p" id="mwAnA">,</span> <span class="n" id="mwAnE">x</span><span class="p" id="mwAnI">:</span> <span class="nb" id="mwAnM">float</span><span class="p" id="mwAnQ">):</span> <span class="bp" id="mwAnU">self</span><span class="o" id="mwAnY">.</span><span class="n" id="mwAnc">count</span> <span class="o" id="mwAng">-=</span> <span class="mi" id="mwAnk">1</span> <span class="n" id="mwAno">new_mean</span> <span class="o" id="mwAns">=</span> <span class="bp" id="mwAnw">self</span><span class="o" id="mwAn0">.</span><span class="n" id="mwAn4">mean</span> <span class="bp" id="mwAn8">self</span><span class="o" id="mwAoA">.</span><span class="n" id="mwAoE">mean</span> <span class="o" id="mwAoI">-=</span> <span class="p" id="mwAoM">(</span><span class="n" id="mwAoQ">x</span> <span class="o" id="mwAoU">-</span> <span class="bp" id="mwAoY">self</span><span class="o" id="mwAoc">.</span><span class="n" id="mwAog">mean</span><span class="p" id="mwAok">)</span> <span class="o" id="mwAoo">/</span> <span class="bp" id="mwAos">self</span><span class="o" id="mwAow">.</span><span class="n" id="mwAo0">count</span> <span class="bp" id="mwAo4">self</span><span class="o" id="mwAo8">.</span><span class="n" id="mwApA">M2</span> <span class="o" id="mwApE">-=</span> <span class="p" id="mwApI">(</span><span class="n" id="mwApM">x</span> <span class="o" id="mwApQ">-</span> <span class="n" id="mwApU">new_mean</span><span class="p" id="mwApY">)</span> <span class="o" id="mwApc">*</span> <span class="p" id="mwApg">(</span><span class="n" id="mwApk">x</span> <span class="o" id="mwApo">-</span> <span class="bp" id="mwAps">self</span><span class="o" id="mwApw">.</span><span class="n" id="mwAp0">mean</span><span class="p" id="mwAp4">)</span> <span class="k" id="mwAp8">def</span><span class="w" id="mwAqA"> </span><span class="nf" id="mwAqE">get_mean</span><span class="p" id="mwAqI">(</span><span class="bp" id="mwAqM">self</span><span class="p" id="mwAqQ">)</span> <span class="o" id="mwAqU">-></span> <span class="nb" id="mwAqY">float</span><span class="p" id="mwAqc">:</span> <span class="k" id="mwAqg">return</span> <span class="bp" id="mwAqk">self</span><span class="o" id="mwAqo">.</span><span class="n" id="mwAqs">mean</span> <span class="k" id="mwAqw">def</span><span class="w" id="mwAq0"> </span><span class="nf" id="mwAq4">get_variance</span><span class="p" id="mwAq8">(</span><span class="bp" id="mwArA">self</span><span class="p" id="mwArE">)</span> <span class="o" id="mwArI">-></span> <span class="nb" id="mwArM">float</span><span class="p" id="mwArQ">:</span> <span class="k" id="mwArU">return</span> <span class="bp" id="mwArY">self</span><span class="o" id="mwArc">.</span><span class="n" id="mwArg">M2</span> <span class="o" id="mwArk">/</span> <span class="bp" id="mwAro">self</span><span class="o" id="mwArs">.</span><span class="n" id="mwArw">count</span> <span class="k" id="mwAr0">def</span><span class="w" id="mwAr4"> </span><span class="nf" id="mwAr8">get_sample_variance</span><span class="p" id="mwAsA">(</span><span class="bp" id="mwAsE">self</span><span class="p" id="mwAsI">)</span> <span class="o" id="mwAsM">-></span> <span class="nb" id="mwAsQ">float</span><span class="p" id="mwAsU">:</span> <span class="k" id="mwAsY">return</span> <span class="bp" id="mwAsc">self</span><span class="o" id="mwAsg">.</span><span class="n" id="mwAsk">M2</span> <span class="o" id="mwAso">/</span> <span class="p" id="mwAss">(</span><span class="bp" id="mwAsw">self</span><span class="o" id="mwAs0">.</span><span class="n" id="mwAs4">count</span> <span class="o" id="mwAs8">-</span> <span class="mi" id="mwAtA">1</span><span class="p" id="mwAtE">)</span>

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) suggests this incremental algorithm:

<span id="mwAuA"></span><span class="kn" id="mwAuE">from</span><span class="w" id="mwAuI"> </span><span class="nn" id="mwAuM">collections</span><span class="w" id="mwAuQ"> </span><span class="kn" id="mwAuU">import</span> <span class="n" id="mwAuY">namedtuple</span> <span class="n" id="mwAuc">WeightedVariances</span> <span class="o" id="mwAug">=</span> <span class="n" id="mwAuk">namedtuple</span><span class="p" id="mwAuo">(</span><span class="s2" id="mwAus">"WeightedVariances"</span><span class="p" id="mwAuw">,</span> <span class="s2" id="mwAu0">"pop freq reli"</span><span class="p" id="mwAu4">)</span> <span class="k" id="mwAu8">def</span><span class="w" id="mwAvA"> </span><span class="nf" id="mwAvE">weighted_incremental_variance</span><span class="p" id="mwAvI">(</span><span class="n" id="mwAvM">data_weight_pairs</span><span class="p" id="mwAvQ">):</span> <span class="n" id="mwAvU">w_sum</span> <span class="o" id="mwAvY">=</span> <span class="n" id="mwAvc">w_sum2</span> <span class="o" id="mwAvg">=</span> <span class="n" id="mwAvk">mean</span> <span class="o" id="mwAvo">=</span> <span class="n" id="mwAvs">S</span> <span class="o" id="mwAvw">=</span> <span class="mi" id="mwAv0">0</span> <span class="k" id="mwAv4">for</span> <span class="n" id="mwAv8">x</span><span class="p" id="mwAwA">,</span> <span class="n" id="mwAwE">w</span> <span class="ow" id="mwAwI">in</span> <span class="n" id="mwAwM">data_weight_pairs</span><span class="p" id="mwAwQ">:</span> <span class="n" id="mwAwU">w_sum</span> <span class="o" id="mwAwY">=</span> <span class="n" id="mwAwc">w_sum</span> <span class="o" id="mwAwg">+</span> <span class="n" id="mwAwk">w</span> <span class="n" id="mwAwo">w_sum2</span> <span class="o" id="mwAws">=</span> <span class="n" id="mwAww">w_sum2</span> <span class="o" id="mwAw0">+</span> <span class="n" id="mwAw4">w</span><span class="o" id="mwAw8">**</span><span class="mi" id="mwAxA">2</span> <span class="n" id="mwAxE">mean_old</span> <span class="o" id="mwAxI">=</span> <span class="n" id="mwAxM">mean</span> <span class="n" id="mwAxQ">mean</span> <span class="o" id="mwAxU">=</span> <span class="n" id="mwAxY">mean_old</span> <span class="o" id="mwAxc">+</span> <span class="p" id="mwAxg">(</span><span class="n" id="mwAxk">w</span> <span class="o" id="mwAxo">/</span> <span class="n" id="mwAxs">w_sum</span><span class="p" id="mwAxw">)</span> <span class="o" id="mwAx0">*</span> <span class="p" id="mwAx4">(</span><span class="n" id="mwAx8">x</span> <span class="o" id="mwAyA">-</span> <span class="n" id="mwAyE">mean_old</span><span class="p" id="mwAyI">)</span> <span class="n" id="mwAyM">S</span> <span class="o" id="mwAyQ">=</span> <span class="n" id="mwAyU">S</span> <span class="o" id="mwAyY">+</span> <span class="n" id="mwAyc">w</span> <span class="o" id="mwAyg">*</span> <span class="p" id="mwAyk">(</span><span class="n" id="mwAyo">x</span> <span class="o" id="mwAys">-</span> <span class="n" id="mwAyw">mean_old</span><span class="p" id="mwAy0">)</span> <span class="o" id="mwAy4">*</span> <span class="p" id="mwAy8">(</span><span class="n" id="mwAzA">x</span> <span class="o" id="mwAzE">-</span> <span class="n" id="mwAzI">mean</span><span class="p" id="mwAzM">)</span> <span class="n" id="mwAzQ">population_variance</span> <span class="o" id="mwAzU">=</span> <span class="n" id="mwAzY">S</span> <span class="o" id="mwAzc">/</span> <span class="n" id="mwAzg">w_sum</span> <span class="c1" id="mwAzk"># Bessel's correction for weighted samples</span> <span class="c1" id="mwAzo"># Frequency weights</span> <span class="n" id="mwAzs">sample_frequency_variance</span> <span class="o" id="mwAzw">=</span> <span class="n" id="mwAz0">S</span> <span class="o" id="mwAz4">/</span> <span class="p" id="mwAz8">(</span><span class="n" id="mwA0A">w_sum</span> <span class="o" id="mwA0E">-</span> <span class="mi" id="mwA0I">1</span><span class="p" id="mwA0M">)</span> <span class="c1" id="mwA0Q"># Reliability weights</span> <span class="n" id="mwA0U">sample_reliability_variance</span> <span class="o" id="mwA0Y">=</span> <span class="n" id="mwA0c">S</span> <span class="o" id="mwA0g">/</span> <span class="p" id="mwA0k">(</span><span class="n" id="mwA0o">w_sum</span> <span class="o" id="mwA0s">-</span> <span class="n" id="mwA0w">w_sum2</span> <span class="o" id="mwA00">/</span> <span class="n" id="mwA04">w_sum</span><span class="p" id="mwA08">)</span> <span class="k" id="mwA1A">return</span> <span class="n" id="mwA1E">WeightedVariances</span><span class="p" id="mwA1I">(</span><span class="n" id="mwA1M">population_variance</span><span class="p" id="mwA1Q">,</span> <span class="n" id="mwA1U">sample_frequency_variance</span><span class="p" id="mwA1Y">,</span> <span class="n" id="mwA1c">sample_reliability_variance</span><span class="p" id="mwA1g">)</span>

Parallel algorithm

Chan et al. note that Welford's online algorithm detailed above is a special case of an algorithm that works for combining arbitrary sets $A$ and $B$ :

\begin{matrix} n_{A B} & = n_{A} + n_{B} \\ δ & = {\bar{x}}_{B} - {\bar{x}}_{A} \\ {\bar{x}}_{A B} & = {\bar{x}}_{A} + δ \cdot \frac{n_{B}}{n_{A B}} \\ M_{2, A B} & = M_{2, A} + M_{2, B} + δ^{2} \cdot \frac{n_{A} n_{B}}{n_{A B}} \end{matrix}

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 $n_{A} \approx n_{B}$ and both are large, because the numerical error in $δ = {\bar{x}}_{B} - {\bar{x}}_{A}$ is not scaled down in the way that it is in the $n_{B} = 1$ case. In such cases, prefer ${\bar{x}}_{A B} = \frac{n_{A} {\bar{x}}_{A} + n_{B} {\bar{x}}_{B}}{n_{A B}}$ . Reusing the WelfordVariance class from above, we have:

<span id="mwA3E"></span><span class="k" id="mwA3I">def</span><span class="w" id="mwA3M"> </span><span class="nf" id="mwA3Q">merge</span><span class="p" id="mwA3U">(</span><span class="n" id="mwA3Y">a</span><span class="p" id="mwA3c">:</span> <span class="n" id="mwA3g">WelfordVariance</span><span class="p" id="mwA3k">,</span> <span class="n" id="mwA3o">b</span><span class="p" id="mwA3s">:</span> <span class="n" id="mwA3w">WelfordVariance</span><span class="p" id="mwA30">)</span> <span class="o" id="mwA34">-></span> <span class="n" id="mwA38">WelfordVariance</span><span class="p" id="mwA4A">:</span> <span class="n" id="mwA4E">ab</span> <span class="o" id="mwA4I">=</span> <span class="n" id="mwA4M">WelfordVariance</span><span class="p" id="mwA4Q">()</span> <span class="n" id="mwA4U">ab</span><span class="o" id="mwA4Y">.</span><span class="n" id="mwA4c">count</span> <span class="o" id="mwA4g">=</span> <span class="n" id="mwA4k">a</span><span class="o" id="mwA4o">.</span><span class="n" id="mwA4s">count</span> <span class="o" id="mwA4w">+</span> <span class="n" id="mwA40">b</span><span class="o" id="mwA44">.</span><span class="n" id="mwA48">count</span> <span class="n" id="mwA5A">delta</span> <span class="o" id="mwA5E">=</span> <span class="n" id="mwA5I">b</span><span class="o" id="mwA5M">.</span><span class="n" id="mwA5Q">mean</span> <span class="o" id="mwA5U">-</span> <span class="n" id="mwA5Y">a</span><span class="o" id="mwA5c">.</span><span class="n" id="mwA5g">mean</span> <span class="n" id="mwA5k">ab</span><span class="o" id="mwA5o">.</span><span class="n" id="mwA5s">mean</span> <span class="o" id="mwA5w">=</span> <span class="p" id="mwA50">(</span><span class="n" id="mwA54">a</span><span class="o" id="mwA58">.</span><span class="n" id="mwA6A">count</span> <span class="o" id="mwA6E">*</span> <span class="n" id="mwA6I">a</span><span class="o" id="mwA6M">.</span><span class="n" id="mwA6Q">mean</span> <span class="o" id="mwA6U">+</span> <span class="n" id="mwA6Y">b</span><span class="o" id="mwA6c">.</span><span class="n" id="mwA6g">count</span> <span class="o" id="mwA6k">*</span> <span class="n" id="mwA6o">b</span><span class="o" id="mwA6s">.</span><span class="n" id="mwA6w">mean</span><span class="p" id="mwA60">)</span> <span class="o" id="mwA64">/</span> <span class="n" id="mwA68">ab</span><span class="o" id="mwA7A">.</span><span class="n" id="mwA7E">count</span> <span class="n" id="mwA7I">ab</span><span class="o" id="mwA7M">.</span><span class="n" id="mwA7Q">M2</span> <span class="o" id="mwA7U">=</span> <span class="n" id="mwA7Y">a</span><span class="o" id="mwA7c">.</span><span class="n" id="mwA7g">M2</span> <span class="o" id="mwA7k">+</span> <span class="n" id="mwA7o">b</span><span class="o" id="mwA7s">.</span><span class="n" id="mwA7w">M2</span> <span class="o" id="mwA70">+</span> <span class="n" id="mwA74">delta</span><span class="o" id="mwA78">**</span><span class="mi" id="mwA8A">2</span> <span class="o" id="mwA8E">*</span> <span class="n" id="mwA8I">a</span><span class="o" id="mwA8M">.</span><span class="n" id="mwA8Q">count</span> <span class="o" id="mwA8U">*</span> <span class="n" id="mwA8Y">b</span><span class="o" id="mwA8c">.</span><span class="n" id="mwA8g">count</span> <span class="o" id="mwA8k">/</span> <span class="n" id="mwA8o">ab</span><span class="o" id="mwA8s">.</span><span class="n" id="mwA8w">count</span> <span class="k" id="mwA80">return</span> <span class="n" id="mwA84">ab</span> <span class="c1" id="mwA88"># example:</span> <span class="n" id="mwA9A">ab</span> <span class="o" id="mwA9E">=</span> <span class="n" id="mwA9I">merge</span><span class="p" id="mwA9M">(</span><span class="n" id="mwA9Q">a</span><span class="p" id="mwA9U">,</span> <span class="n" id="mwA9Y">b</span><span class="p" id="mwA9c">)</span> <span class="nb" id="mwA9g">print</span><span class="p" id="mwA9k">(</span><span class="n" id="mwA9o">ab</span><span class="o" id="mwA9s">.</span><span class="n" id="mwA9w">get_sample_variance</span><span class="p" id="mwA90">())</span>

This algorithm allows one to divide a dataset into several pieces, run them in parallel, and then merge the results together. This enables parallelization in any way, including AVX, with GPUs, and computer clusters. The algorithm can also be adapted for higher-order statistics as well as 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 (10⁸ + 4, 10⁸ + 7, 10⁸ + 13, 10⁸ + 16), 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 (10⁹ + 4, 10⁹ + 7, 10⁹ + 13, 10⁹ + 16). 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:

\begin{matrix} M_{3, X} = M_{3, A} + M_{3, B} & + δ^{3} \frac{n_{A} n_{B} (n_{A} - n_{B})}{n_{X}^{2}} + 3 δ \frac{n_{A} M_{2, B} - n_{B} M_{2, A}}{n_{X}} \\ M_{4, X} = M_{4, A} + M_{4, B} & + δ^{4} \frac{n_{A} n_{B} (n_{A}^{2} - n_{A} n_{B} + n_{B}^{2})}{n_{X}^{3}} \\ + 6 δ^{2} \frac{n_{A}^{2} M_{2, B} + n_{B}^{2} M_{2, A}}{n_{X}^{2}} + 4 δ \frac{n_{A} M_{3, B} - n_{B} M_{3, A}}{n_{X}} \end{matrix}

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

\begin{matrix} skewness = g_{1} = \frac{\sqrt{n} M_{3}}{M_{2}^{3 / 2}}, \\ kurtosis = g_{2} = \frac{n M_{4}}{M_{2}^{2}} - 3. \end{matrix}

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

\begin{matrix} δ & = x - m \\ m^{'} & = m + \frac{δ}{n} \\ M_{2}^{'} & = M_{2} + δ^{2} \frac{n - 1}{n} \\ M_{3}^{'} & = M_{3} + δ^{3} \frac{(n - 1) (n - 2)}{n^{2}} - \frac{3 δ M_{2}}{n} \\ M_{4}^{'} & = M_{4} + \frac{δ^{4} (n - 1) (n^{2} - 3 n + 3)}{n^{3}} + \frac{6 δ^{2} M_{2}}{n^{2}} - \frac{4 δ M_{3}}{n} \end{matrix}

By preserving the value $δ / 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:

<span id="mwBBU"></span><span class="k" id="mwBBY">def</span><span class="w" id="mwBBc"> </span><span class="nf" id="mwBBg">online_kurtosis</span><span class="p" id="mwBBk">(</span><span class="n" id="mwBBo">data</span><span class="p" id="mwBBs">):</span> <span class="n" id="mwBBw">n</span> <span class="o" id="mwBB0">=</span> <span class="n" id="mwBB4">mean</span> <span class="o" id="mwBB8">=</span> <span class="n" id="mwBCA">M2</span> <span class="o" id="mwBCE">=</span> <span class="n" id="mwBCI">M3</span> <span class="o" id="mwBCM">=</span> <span class="n" id="mwBCQ">M4</span> <span class="o" id="mwBCU">=</span> <span class="mi" id="mwBCY">0</span> <span class="k" id="mwBCc">for</span> <span class="n" id="mwBCg">x</span> <span class="ow" id="mwBCk">in</span> <span class="n" id="mwBCo">data</span><span class="p" id="mwBCs">:</span> <span class="n" id="mwBCw">n1</span> <span class="o" id="mwBC0">=</span> <span class="n" id="mwBC4">n</span> <span class="n" id="mwBC8">n</span> <span class="o" id="mwBDA">=</span> <span class="n" id="mwBDE">n</span> <span class="o" id="mwBDI">+</span> <span class="mi" id="mwBDM">1</span> <span class="n" id="mwBDQ">delta</span> <span class="o" id="mwBDU">=</span> <span class="n" id="mwBDY">x</span> <span class="o" id="mwBDc">-</span> <span class="n" id="mwBDg">mean</span> <span class="n" id="mwBDk">delta_n</span> <span class="o" id="mwBDo">=</span> <span class="n" id="mwBDs">delta</span> <span class="o" id="mwBDw">/</span> <span class="n" id="mwBD0">n</span> <span class="n" id="mwBD4">delta_n2</span> <span class="o" id="mwBD8">=</span> <span class="n" id="mwBEA">delta_n</span><span class="o" id="mwBEE">**</span><span class="mi" id="mwBEI">2</span> <span class="n" id="mwBEM">term1</span> <span class="o" id="mwBEQ">=</span> <span class="n" id="mwBEU">delta</span> <span class="o" id="mwBEY">*</span> <span class="n" id="mwBEc">delta_n</span> <span class="o" id="mwBEg">*</span> <span class="n" id="mwBEk">n1</span> <span class="n" id="mwBEo">mean</span> <span class="o" id="mwBEs">=</span> <span class="n" id="mwBEw">mean</span> <span class="o" id="mwBE0">+</span> <span class="n" id="mwBE4">delta_n</span> <span class="n" id="mwBE8">M4</span> <span class="o" id="mwBFA">=</span> <span class="n" id="mwBFE">M4</span> <span class="o" id="mwBFI">+</span> <span class="n" id="mwBFM">term1</span> <span class="o" id="mwBFQ">*</span> <span class="n" id="mwBFU">delta_n2</span> <span class="o" id="mwBFY">*</span> <span class="p" id="mwBFc">(</span><span class="n" id="mwBFg">n</span><span class="o" id="mwBFk">**</span><span class="mi" id="mwBFo">2</span> <span class="o" id="mwBFs">-</span> <span class="mi" id="mwBFw">3</span><span class="o" id="mwBF0">*</span><span class="n" id="mwBF4">n</span> <span class="o" id="mwBF8">+</span> <span class="mi" id="mwBGA">3</span><span class="p" id="mwBGE">)</span> <span class="o" id="mwBGI">+</span> <span class="mi" id="mwBGM">6</span> <span class="o" id="mwBGQ">*</span> <span class="n" id="mwBGU">delta_n2</span> <span class="o" id="mwBGY">*</span> <span class="n" id="mwBGc">M2</span> <span class="o" id="mwBGg">-</span> <span class="mi" id="mwBGk">4</span> <span class="o" id="mwBGo">*</span> <span class="n" id="mwBGs">delta_n</span> <span class="o" id="mwBGw">*</span> <span class="n" id="mwBG0">M3</span> <span class="n" id="mwBG4">M3</span> <span class="o" id="mwBG8">=</span> <span class="n" id="mwBHA">M3</span> <span class="o" id="mwBHE">+</span> <span class="n" id="mwBHI">term1</span> <span class="o" id="mwBHM">*</span> <span class="n" id="mwBHQ">delta_n</span> <span class="o" id="mwBHU">*</span> <span class="p" id="mwBHY">(</span><span class="n" id="mwBHc">n</span> <span class="o" id="mwBHg">-</span> <span class="mi" id="mwBHk">2</span><span class="p" id="mwBHo">)</span> <span class="o" id="mwBHs">-</span> <span class="mi" id="mwBHw">3</span> <span class="o" id="mwBH0">*</span> <span class="n" id="mwBH4">delta_n</span> <span class="o" id="mwBH8">*</span> <span class="n" id="mwBIA">M2</span> <span class="n" id="mwBIE">M2</span> <span class="o" id="mwBII">=</span> <span class="n" id="mwBIM">M2</span> <span class="o" id="mwBIQ">+</span> <span class="n" id="mwBIU">term1</span> <span class="c1" id="mwBIY"># Note, you may also calculate variance using M2, and skewness using M3</span> <span class="c1" id="mwBIc"># Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.</span> <span class="n" id="mwBIg">kurtosis</span> <span class="o" id="mwBIk">=</span> <span class="p" id="mwBIo">(</span><span class="n" id="mwBIs">n</span> <span class="o" id="mwBIw">*</span> <span class="n" id="mwBI0">M4</span><span class="p" id="mwBI4">)</span> <span class="o" id="mwBI8">/</span> <span class="p" id="mwBJA">(</span><span class="n" id="mwBJE">M2</span><span class="o" id="mwBJI">**</span><span class="mi" id="mwBJM">2</span><span class="p" id="mwBJQ">)</span> <span class="o" id="mwBJU">-</span> <span class="mi" id="mwBJY">3</span> <span class="k" id="mwBJc">return</span> <span class="n" id="mwBJg">kurtosis</span>

Pébaÿ further extends these results to arbitrary-order central moments, for the incremental and the pairwise cases, and subsequently Pébaÿ et al. for weighted and compound moments. One can also find there similar formulas for covariance.

Choi and Sweetman 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:

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

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

m_{n}^{(h)} = \sum_{k = 1}^{K} x_{k}^{n} H (x_{k}) Δ x_{k} = \frac{1}{A} \sum_{k = 1}^{K} x_{k}^{n} h (x_{k}) Δ x_{k}

θ_{n}^{(h)} = \sum_{k = 1}^{K} (x_{k} - m_{1}^{(h)})^{n} H (x_{k}) Δ x_{k} = \frac{1}{A} \sum_{k = 1}^{K} (x_{k} - m_{1}^{(h)})^{n} h (x_{k}) Δ x_{k}

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

m_{n}^{(h)} = \frac{1}{I} \sum_{k = 1}^{K} x_{k}^{n} h (x_{k})

θ_{n}^{(h)} = \frac{1}{I} \sum_{k = 1}^{K} (x_{k} - m_{1}^{(h)})^{n} h (x_{k})

The second approach from Choi and Sweetman 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 $Q$ sets of statistical moments are known: $(γ_{0, q}, μ_{q}, σ_{q}^{2}, α_{3, q}, α_{4, q})$ for $q = 1, 2, \dots, Q$ , then each $γ_{n}$ can be expressed in terms of the equivalent $n$ raw moments:

γ_{n, q} = m_{n, q} γ_{0, q} for n = 1, 2, 3, 4 and q = 1, 2, \dots, Q

where $γ_{0, q}$ is generally taken to be the duration of the $q^{t h}$ time-history, or the number of points if $Δ t$ is constant.

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

γ_{n, c} = \sum_{q = 1}^{Q} γ_{n, q} for n = 0, 1, 2, 3, 4

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

m_{n, c} = \frac{γ_{n, c}}{γ_{0, c}} for n = 1, 2, 3, 4

Known relationships between the raw moments ( $m_{n}$ ) and the central moments ( $θ_{n} = E [(x - μ)^{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:

μ_{c} = m_{1, c} σ_{c}^{2} = θ_{2, c} α_{3, c} = \frac{θ_{3, c}}{σ_{c}^{3}} α_{4, c} = \frac{θ_{4, c}}{σ_{c}^{4}} - 3

Covariance

Very similar algorithms can be used to compute the covariance.

Naïve algorithm

The naïve algorithm is

Cov (X, Y) = \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:

<span id="mwBPg"></span><span class="k" id="mwBPk">def</span><span class="w" id="mwBPo"> </span><span class="nf" id="mwBPs">naive_covariance</span><span class="p" id="mwBPw">(</span><span class="n" id="mwBP0">data1</span><span class="p" id="mwBP4">,</span> <span class="n" id="mwBP8">data2</span><span class="p" id="mwBQA">):</span> <span class="n" id="mwBQE">n</span> <span class="o" id="mwBQI">=</span> <span class="nb" id="mwBQM">len</span><span class="p" id="mwBQQ">(</span><span class="n" id="mwBQU">data1</span><span class="p" id="mwBQY">)</span> <span class="n" id="mwBQc">sum1</span> <span class="o" id="mwBQg">=</span> <span class="nb" id="mwBQk">sum</span><span class="p" id="mwBQo">(</span><span class="n" id="mwBQs">data1</span><span class="p" id="mwBQw">)</span> <span class="n" id="mwBQ0">sum2</span> <span class="o" id="mwBQ4">=</span> <span class="nb" id="mwBQ8">sum</span><span class="p" id="mwBRA">(</span><span class="n" id="mwBRE">data2</span><span class="p" id="mwBRI">)</span> <span class="n" id="mwBRM">sum12</span> <span class="o" id="mwBRQ">=</span> <span class="nb" id="mwBRU">sum</span><span class="p" id="mwBRY">([</span><span class="n" id="mwBRc">i1</span> <span class="o" id="mwBRg">*</span> <span class="n" id="mwBRk">i2</span> <span class="k" id="mwBRo">for</span> <span class="n" id="mwBRs">i1</span><span class="p" id="mwBRw">,</span> <span class="n" id="mwBR0">i2</span> <span class="ow" id="mwBR4">in</span> <span class="nb" id="mwBR8">zip</span><span class="p" id="mwBSA">(</span><span class="n" id="mwBSE">data1</span><span class="p" id="mwBSI">,</span> <span class="n" id="mwBSM">data2</span><span class="p" id="mwBSQ">)])</span> <span class="n" id="mwBSU">covariance</span> <span class="o" id="mwBSY">=</span> <span class="p" id="mwBSc">(</span><span class="n" id="mwBSg">sum12</span> <span class="o" id="mwBSk">-</span> <span class="n" id="mwBSo">sum1</span> <span class="o" id="mwBSs">*</span> <span class="n" id="mwBSw">sum2</span> <span class="o" id="mwBS0">/</span> <span class="n" id="mwBS4">n</span><span class="p" id="mwBS8">)</span> <span class="o" id="mwBTA">/</span> <span class="n" id="mwBTE">n</span> <span class="k" id="mwBTI">return</span> <span class="n" id="mwBTM">covariance</span>

With estimate of the mean

As for the variance, the covariance of two random variables is also shift-invariant, so given any two constant values $k_{x}$ and $k_{y},$ it can be written:

Cov (X, Y) = Cov (X - k_{x}, Y - k_{y}) = \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:

<span id="mwBT8"></span><span class="k" id="mwBUA">def</span><span class="w" id="mwBUE"> </span><span class="nf" id="mwBUI">shifted_data_covariance</span><span class="p" id="mwBUM">(</span><span class="n" id="mwBUQ">data_x</span><span class="p" id="mwBUU">,</span> <span class="n" id="mwBUY">data_y</span><span class="p" id="mwBUc">):</span> <span class="n" id="mwBUg">n</span> <span class="o" id="mwBUk">=</span> <span class="nb" id="mwBUo">len</span><span class="p" id="mwBUs">(</span><span class="n" id="mwBUw">data_x</span><span class="p" id="mwBU0">)</span> <span class="k" id="mwBU4">if</span> <span class="n" id="mwBU8">n</span> <span class="o" id="mwBVA"><</span> <span class="mi" id="mwBVE">2</span><span class="p" id="mwBVI">:</span> <span class="k" id="mwBVM">return</span> <span class="mi" id="mwBVQ">0</span> <span class="n" id="mwBVU">kx</span> <span class="o" id="mwBVY">=</span> <span class="n" id="mwBVc">data_x</span><span class="p" id="mwBVg">[</span><span class="mi" id="mwBVk">0</span><span class="p" id="mwBVo">]</span> <span class="n" id="mwBVs">ky</span> <span class="o" id="mwBVw">=</span> <span class="n" id="mwBV0">data_y</span><span class="p" id="mwBV4">[</span><span class="mi" id="mwBV8">0</span><span class="p" id="mwBWA">]</span> <span class="n" id="mwBWE">Ex</span> <span class="o" id="mwBWI">=</span> <span class="n" id="mwBWM">Ey</span> <span class="o" id="mwBWQ">=</span> <span class="n" id="mwBWU">Exy</span> <span class="o" id="mwBWY">=</span> <span class="mi" id="mwBWc">0</span> <span class="k" id="mwBWg">for</span> <span class="n" id="mwBWk">ix</span><span class="p" id="mwBWo">,</span> <span class="n" id="mwBWs">iy</span> <span class="ow" id="mwBWw">in</span> <span class="nb" id="mwBW0">zip</span><span class="p" id="mwBW4">(</span><span class="n" id="mwBW8">data_x</span><span class="p" id="mwBXA">,</span> <span class="n" id="mwBXE">data_y</span><span class="p" id="mwBXI">):</span> <span class="c1" id="mwBXM"># for higher precision, use sum():</span> <span class="n" id="mwBXQ">Ex</span> <span class="o" id="mwBXU">+=</span> <span class="n" id="mwBXY">ix</span> <span class="o" id="mwBXc">-</span> <span class="n" id="mwBXg">kx</span> <span class="c1" id="mwBXk"># Ex = sum(ix - kx for ix in data_x) # or sum(ix) - kx * n</span> <span class="n" id="mwBXo">Ey</span> <span class="o" id="mwBXs">+=</span> <span class="n" id="mwBXw">iy</span> <span class="o" id="mwBX0">-</span> <span class="n" id="mwBX4">ky</span> <span class="c1" id="mwBX8"># Ey = sum(iy - ky for iy in data_y) # or sum(iy) - ky * n</span> <span class="n" id="mwBYA">Exy</span> <span class="o" id="mwBYE">+=</span> <span class="p" id="mwBYI">(</span><span class="n" id="mwBYM">ix</span> <span class="o" id="mwBYQ">-</span> <span class="n" id="mwBYU">kx</span><span class="p" id="mwBYY">)</span> <span class="o" id="mwBYc">*</span> <span class="p" id="mwBYg">(</span><span class="n" id="mwBYk">iy</span> <span class="o" id="mwBYo">-</span> <span class="n" id="mwBYs">ky</span><span class="p" id="mwBYw">)</span> <span class="c1" id="mwBY0"># Exy = sum((ix - kx) * (iy - ky) for ix, iy in zip(data_x, data_y)) </span> <span class="k" id="mwBY4">return</span> <span class="p" id="mwBY8">(</span><span class="n" id="mwBZA">Exy</span> <span class="o" id="mwBZE">-</span> <span class="n" id="mwBZI">Ex</span> <span class="o" id="mwBZM">*</span> <span class="n" id="mwBZQ">Ey</span> <span class="o" id="mwBZU">/</span> <span class="n" id="mwBZY">n</span><span class="p" id="mwBZc">)</span> <span class="o" id="mwBZg">/</span> <span class="n" id="mwBZk">n</span>

Two-pass

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

\bar{x} = \sum_{i = 1}^{n} x_{i} / n

\bar{y} = \sum_{i = 1}^{n} y_{i} / n

Cov (X, Y) = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y})}{n} .

The two-pass algorithm may be written as:

<span id="mwBaY"></span><span class="k" id="mwBac">def</span><span class="w" id="mwBag"> </span><span class="nf" id="mwBak">two_pass_covariance</span><span class="p" id="mwBao">(</span><span class="n" id="mwBas">data1</span><span class="p" id="mwBaw">,</span> <span class="n" id="mwBa0">data2</span><span class="p" id="mwBa4">):</span> <span class="n" id="mwBa8">n</span> <span class="o" id="mwBbA">=</span> <span class="nb" id="mwBbE">len</span><span class="p" id="mwBbI">(</span><span class="n" id="mwBbM">data1</span><span class="p" id="mwBbQ">)</span> <span class="n" id="mwBbU">mean1</span> <span class="o" id="mwBbY">=</span> <span class="nb" id="mwBbc">sum</span><span class="p" id="mwBbg">(</span><span class="n" id="mwBbk">data1</span><span class="p" id="mwBbo">)</span> <span class="o" id="mwBbs">/</span> <span class="n" id="mwBbw">n</span> <span class="n" id="mwBb0">mean2</span> <span class="o" id="mwBb4">=</span> <span class="nb" id="mwBb8">sum</span><span class="p" id="mwBcA">(</span><span class="n" id="mwBcE">data2</span><span class="p" id="mwBcI">)</span> <span class="o" id="mwBcM">/</span> <span class="n" id="mwBcQ">n</span> <span class="n" id="mwBcU">covariance</span> <span class="o" id="mwBcY">=</span> <span class="mi" id="mwBcc">0</span> <span class="k" id="mwBcg">for</span> <span class="n" id="mwBck">i1</span><span class="p" id="mwBco">,</span> <span class="n" id="mwBcs">i2</span> <span class="ow" id="mwBcw">in</span> <span class="nb" id="mwBc0">zip</span><span class="p" id="mwBc4">(</span><span class="n" id="mwBc8">data1</span><span class="p" id="mwBdA">,</span> <span class="n" id="mwBdE">data2</span><span class="p" id="mwBdI">):</span> <span class="c1" id="mwBdM"># for higher precision, use sum():</span> <span class="n" id="mwBdQ">a</span> <span class="o" id="mwBdU">=</span> <span class="n" id="mwBdY">i1</span> <span class="o" id="mwBdc">-</span> <span class="n" id="mwBdg">mean1</span> <span class="c1" id="mwBdk"># covariance = sum((i1 - mean1) * (i2 - mean2) for i1, i2 in zip(data1, data2))</span> <span class="n" id="mwBdo">b</span> <span class="o" id="mwBds">=</span> <span class="n" id="mwBdw">i2</span> <span class="o" id="mwBd0">-</span> <span class="n" id="mwBd4">mean2</span> <span class="n" id="mwBd8">covariance</span> <span class="o" id="mwBeA">+=</span> <span class="n" id="mwBeE">a</span> <span class="o" id="mwBeI">*</span> <span class="n" id="mwBeM">b</span> <span class="o" id="mwBeQ">/</span> <span class="n" id="mwBeU">n</span> <span class="k" id="mwBeY">return</span> <span class="n" id="mwBec">covariance</span>

Online

A stable one-pass algorithm exists, similar to the online algorithm for computing the variance, that computes co-moment $C_{n} = \sum_{i = 1}^{n} (x_{i} - {\bar{x}}_{n}) (y_{i} - {\bar{y}}_{n})$ :

\begin{aligned} {\bar{x}}_{n} & = {\bar{x}}_{n - 1} & + & \frac{x_{n} - {\bar{x}}_{n - 1}}{n} \\ {\bar{y}}_{n} & = {\bar{y}}_{n - 1} & + & \frac{y_{n} - {\bar{y}}_{n - 1}}{n} \\ C_{n} & = C_{n - 1} & + & (x_{n} - {\bar{x}}_{n}) (y_{n} - {\bar{y}}_{n - 1}) \\ = C_{n - 1} & + & (x_{n} - {\bar{x}}_{n - 1}) (y_{n} - {\bar{y}}_{n}) \end{aligned}

The apparent asymmetry in that last equation is due to the fact that $(x_{n} - {\bar{x}}_{n}) = \frac{n - 1}{n} (x_{n} - {\bar{x}}_{n - 1})$ , so both update terms are equal to $\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 the covariance can be computed as

\begin{matrix} {Cov}_{N} (X, Y) = \frac{C_{N}}{N} & = \frac{{Cov}_{N - 1} (X, Y) \cdot (N - 1) + (x_{n} - {\bar{x}}_{n}) (y_{n} - {\bar{y}}_{n - 1})}{N} \\ = \frac{{Cov}_{N - 1} (X, Y) \cdot (N - 1) + (x_{n} - {\bar{x}}_{n - 1}) (y_{n} - {\bar{y}}_{n})}{N} \\ = \frac{{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} \\ = \frac{{Cov}_{N - 1} (X, Y) \cdot (N - 1) + \frac{N}{N - 1} (x_{n} - {\bar{x}}_{n}) (y_{n} - {\bar{y}}_{n})}{N} . \end{matrix}

<span id="mwBfc"></span><span class="k" id="mwBfg">def</span><span class="w" id="mwBfk"> </span><span class="nf" id="mwBfo">online_covariance</span><span class="p" id="mwBfs">(</span><span class="n" id="mwBfw">data1</span><span class="p" id="mwBf0">,</span> <span class="n" id="mwBf4">data2</span><span class="p" id="mwBf8">):</span> <span class="n" id="mwBgA">meanx</span> <span class="o" id="mwBgE">=</span> <span class="n" id="mwBgI">meany</span> <span class="o" id="mwBgM">=</span> <span class="n" id="mwBgQ">C</span> <span class="o" id="mwBgU">=</span> <span class="n" id="mwBgY">n</span> <span class="o" id="mwBgc">=</span> <span class="mi" id="mwBgg">0</span> <span class="k" id="mwBgk">for</span> <span class="n" id="mwBgo">x</span><span class="p" id="mwBgs">,</span> <span class="n" id="mwBgw">y</span> <span class="ow" id="mwBg0">in</span> <span class="nb" id="mwBg4">zip</span><span class="p" id="mwBg8">(</span><span class="n" id="mwBhA">data1</span><span class="p" id="mwBhE">,</span> <span class="n" id="mwBhI">data2</span><span class="p" id="mwBhM">):</span> <span class="n" id="mwBhQ">n</span> <span class="o" id="mwBhU">+=</span> <span class="mi" id="mwBhY">1</span> <span class="n" id="mwBhc">dx</span> <span class="o" id="mwBhg">=</span> <span class="n" id="mwBhk">x</span> <span class="o" id="mwBho">-</span> <span class="n" id="mwBhs">meanx</span> <span class="n" id="mwBhw">meanx</span> <span class="o" id="mwBh0">+=</span> <span class="n" id="mwBh4">dx</span> <span class="o" id="mwBh8">/</span> <span class="n" id="mwBiA">n</span> <span class="n" id="mwBiE">meany</span> <span class="o" id="mwBiI">+=</span> <span class="p" id="mwBiM">(</span><span class="n" id="mwBiQ">y</span> <span class="o" id="mwBiU">-</span> <span class="n" id="mwBiY">meany</span><span class="p" id="mwBic">)</span> <span class="o" id="mwBig">/</span> <span class="n" id="mwBik">n</span> <span class="n" id="mwBio">C</span> <span class="o" id="mwBis">+=</span> <span class="n" id="mwBiw">dx</span> <span class="o" id="mwBi0">*</span> <span class="p" id="mwBi4">(</span><span class="n" id="mwBi8">y</span> <span class="o" id="mwBjA">-</span> <span class="n" id="mwBjE">meany</span><span class="p" id="mwBjI">)</span> <span class="c1" id="mwBjM"># covariance and Bessel's correction for sample</span> <span class="k" id="mwBjQ">return</span> <span class="p" id="mwBjU">(</span><span class="n" id="mwBjY">C</span> <span class="o" id="mwBjc">/</span> <span class="n" id="mwBjg">n</span><span class="p" id="mwBjk">,</span> <span class="n" id="mwBjo">C</span> <span class="o" id="mwBjs">/</span> <span class="p" id="mwBjw">(</span><span class="n" id="mwBj0">n</span> <span class="o" id="mwBj4">-</span> <span class="mi" id="mwBj8">1</span><span class="p" id="mwBkA">))</span>

A small modification can also be made to compute the weighted covariance:

<span id="mwBkQ"></span><span class="k" id="mwBkU">def</span><span class="w" id="mwBkY"> </span><span class="nf" id="mwBkc">online_weighted_covariance</span><span class="p" id="mwBkg">(</span><span class="n" id="mwBkk">data1</span><span class="p" id="mwBko">,</span> <span class="n" id="mwBks">data2</span><span class="p" id="mwBkw">,</span> <span class="n" id="mwBk0">data3</span><span class="p" id="mwBk4">):</span> <span class="n" id="mwBk8">meanx</span> <span class="o" id="mwBlA">=</span> <span class="n" id="mwBlE">meany</span> <span class="o" id="mwBlI">=</span> <span class="mi" id="mwBlM">0</span> <span class="n" id="mwBlQ">wsum</span> <span class="o" id="mwBlU">=</span> <span class="n" id="mwBlY">wsum2</span> <span class="o" id="mwBlc">=</span> <span class="mi" id="mwBlg">0</span> <span class="n" id="mwBlk">C</span> <span class="o" id="mwBlo">=</span> <span class="mi" id="mwBls">0</span> <span class="k" id="mwBlw">for</span> <span class="n" id="mwBl0">x</span><span class="p" id="mwBl4">,</span> <span class="n" id="mwBl8">y</span><span class="p" id="mwBmA">,</span> <span class="n" id="mwBmE">w</span> <span class="ow" id="mwBmI">in</span> <span class="nb" id="mwBmM">zip</span><span class="p" id="mwBmQ">(</span><span class="n" id="mwBmU">data1</span><span class="p" id="mwBmY">,</span> <span class="n" id="mwBmc">data2</span><span class="p" id="mwBmg">,</span> <span class="n" id="mwBmk">data3</span><span class="p" id="mwBmo">):</span> <span class="n" id="mwBms">wsum</span> <span class="o" id="mwBmw">+=</span> <span class="n" id="mwBm0">w</span> <span class="n" id="mwBm4">wsum2</span> <span class="o" id="mwBm8">+=</span> <span class="n" id="mwBnA">w</span> <span class="o" id="mwBnE">*</span> <span class="n" id="mwBnI">w</span> <span class="n" id="mwBnM">dx</span> <span class="o" id="mwBnQ">=</span> <span class="n" id="mwBnU">x</span> <span class="o" id="mwBnY">-</span> <span class="n" id="mwBnc">meanx</span> <span class="n" id="mwBng">meanx</span> <span class="o" id="mwBnk">+=</span> <span class="p" id="mwBno">(</span><span class="n" id="mwBns">w</span> <span class="o" id="mwBnw">/</span> <span class="n" id="mwBn0">wsum</span><span class="p" id="mwBn4">)</span> <span class="o" id="mwBn8">*</span> <span class="n" id="mwBoA">dx</span> <span class="n" id="mwBoE">meany</span> <span class="o" id="mwBoI">+=</span> <span class="p" id="mwBoM">(</span><span class="n" id="mwBoQ">w</span> <span class="o" id="mwBoU">/</span> <span class="n" id="mwBoY">wsum</span><span class="p" id="mwBoc">)</span> <span class="o" id="mwBog">*</span> <span class="p" id="mwBok">(</span><span class="n" id="mwBoo">y</span> <span class="o" id="mwBos">-</span> <span class="n" id="mwBow">meany</span><span class="p" id="mwBo0">)</span> <span class="n" id="mwBo4">C</span> <span class="o" id="mwBo8">+=</span> <span class="n" id="mwBpA">w</span> <span class="o" id="mwBpE">*</span> <span class="n" id="mwBpI">dx</span> <span class="o" id="mwBpM">*</span> <span class="p" id="mwBpQ">(</span><span class="n" id="mwBpU">y</span> <span class="o" id="mwBpY">-</span> <span class="n" id="mwBpc">meany</span><span class="p" id="mwBpg">)</span> <span class="n" id="mwBpk">population_covar</span> <span class="o" id="mwBpo">=</span> <span class="n" id="mwBps">C</span> <span class="o" id="mwBpw">/</span> <span class="n" id="mwBp0">wsum</span> <span class="c1" id="mwBp4"># Bessel's correction for sample variance</span> <span class="c1" id="mwBp8"># Frequency weights</span> <span class="n" id="mwBqA">sample_frequency_covar</span> <span class="o" id="mwBqE">=</span> <span class="n" id="mwBqI">C</span> <span class="o" id="mwBqM">/</span> <span class="p" id="mwBqQ">(</span><span class="n" id="mwBqU">wsum</span> <span class="o" id="mwBqY">-</span> <span class="mi" id="mwBqc">1</span><span class="p" id="mwBqg">)</span> <span class="c1" id="mwBqk"># Reliability weights</span> <span class="n" id="mwBqo">sample_reliability_covar</span> <span class="o" id="mwBqs">=</span> <span class="n" id="mwBqw">C</span> <span class="o" id="mwBq0">/</span> <span class="p" id="mwBq4">(</span><span class="n" id="mwBq8">wsum</span> <span class="o" id="mwBrA">-</span> <span class="n" id="mwBrE">wsum2</span> <span class="o" id="mwBrI">/</span> <span class="n" id="mwBrM">wsum</span><span class="p" id="mwBrQ">)</span>

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

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}} .

Weighted batched version

A version of the weighted online algorithm that does batched updated also exists: let $w_{1}, \dots w_{N}$ denote the weights, and write

\begin{aligned} {\bar{x}}_{n + k} & = {\bar{x}}_{n} & + & \frac{\sum_{i = n + 1}^{n + k} w_{i} (x_{i} - {\bar{x}}_{n})}{\sum_{i = 1}^{n + k} w_{i}} \\ {\bar{y}}_{n + k} & = {\bar{y}}_{n} & + & \frac{\sum_{i = n + 1}^{n + k} w_{i} (y_{i} - {\bar{y}}_{n})}{\sum_{i = 1}^{n + k} w_{i}} \\ C_{n + k} & = C_{n} & + & \sum_{i = n + 1}^{n + k} w_{i} (x_{i} - {\bar{x}}_{n + k}) (y_{i} - {\bar{y}}_{n}) \\ = C_{n} & + & \sum_{i = n + 1}^{n + k} w_{i} (x_{i} - {\bar{x}}_{n}) (y_{i} - {\bar{y}}_{n + k}) \end{aligned}

The covariance can then be computed as

{Cov}_{N} (X, Y) = \frac{C_{N}}{\sum_{i = 1}^{N} w_{i}}

(This can be used with Python's sum for better accuracy. Block updating is also related to the parallel/merge algorithm.)

Naïve algorithm

Computing shifted data

Two-pass algorithm

Incremental algorithms

Incremental shifted data

Welford's online algorithm

Weighted incremental algorithm

Parallel algorithm

Example

Higher-order statistics

Covariance

Naïve algorithm

With estimate of the mean

Two-pass

Online

Weighted batched version

More in Topics in mathematics