Welford's Algorithm: Running Variance in Python

Values arrive one at a time in xs. After each value, record [mean, variance] of everything seen so far, where the variance is the sample variance (divide by n − 1), defined as 0.0 when n = 1. Return the list of records.

Use O(1) memory beyond the output: you may not keep the values. The answer must be accurate for values of size 10⁹ that differ by single units.

Examples

Input xs = [2, 4, 6]

Output [[2, 0], [3, 2], [4, 4]]

After 2: (2, 0). After 4: (3, 2). After 6: (4, 4).

Input xs = [1.000e9, 1.000e9, 1.000e9, 1.000e9]

Output [[1.000e9, 0], [1.000e9, 4.5], [1.000e9, 21], [1.000e9, 30]]

The final variance is 30. E[x²] − E[x]² in float64 gets this badly wrong.

Input xs = [5]

Output [[5, 0]]

A single value has variance 0 by definition here.

Constraints

1 ≤ len(xs) ≤ 10⁵
|xs[i]| ≤ 10¹²
Checked to a relative tolerance of 10⁻⁶

Hints

Hint 1

Keep three numbers: the count, the mean, and the sum of squared deviations from the mean.

Hint 2

When x arrives, the mean moves by (x − mean)/n. How does the sum of squared deviations change?

Hint 3

m2 += (x − old_mean) · (x − new_mean). Then the variance is m2/(n − 1).

Where candidates go wrong

Keeping Σx and Σx² and computing Σx²/n − mean². Exact in algebra, wrong by hundreds for the second example.
Dividing by n instead of n − 1.
Updating m2 with delta² (the old deviation squared) instead of the product of old and new deviations.

The idea

Let δ = x − mean_old. The new mean is mean_old + δ/n. The sum of squared deviations about the mean grows by exactly δ·(x − mean_new), which you can check by expanding both sums. Every quantity in the update is a deviation, on the scale of the spread rather than the size of the values, so nothing cancels catastrophically. The error is a few units in the last place of the deviations, not of x².

Time O(n) · Space O(1) beyond the output

Reference solution

def running_stats(xs):
    out = []
    n, mean, m2 = 0, 0.0, 0.0           # m2: sum of squared deviations from the current mean
    for x in xs:
        n += 1
        delta = x - mean
        mean += delta / n
        m2 += delta * (x - mean)        # old deviation times new deviation
        out.append([mean, m2 / (n - 1) if n > 1 else 0.0])
    return out

Console⌘/Ctrl + Enter runs

Python runs in your browser with NumPy. The first run downloads the runtime, which takes a few seconds.

Why interviewers ask it

Trading firms ask this one often, because streaming statistics run all through their systems and the naive formula breaks on exactly the data they handle: large prices that move by small amounts.

The theory behind it is in the lesson Statistics coding questions: variance, AUC and rank correlation.

Keep going

Easy12 min

Running mean and variance in one pass