Online Mean and Variance with Welford's Algorithm

1. From Formulas to an Algorithm

Let's start by recalling the two core recurrence relations we derived. For a new data point $x_k$, the mean ($\mu$) and the sum of squared differences from the mean ($M$) are updated as follows:

1. Mean update: ${\mu }_k={\mu }_{k-1}+\frac{x_k-{\mu }_{k-1}}{k}$
2. Sum of squares update: $M_k=M_{k-1}+(x_k-{\mu }_{k-1})(x_k-{\mu }_k)$

To implement this, we need to maintain state as we process data one point at a time. An algorithm for processing a new value x would look like this:

• State required:

  • count (the number of data points processed so far, $k$)
  • mean (the running mean, ${\mu }_k$)
  • M2 (the running sum of squared differences, $M_k$)

• Initialization: Before processing any data, count, mean, and M2 are all initialized to 0.

• Update Step (for each new value x):

  1. Increment count.
  2. Store the current mean as old_mean.
  3. Calculate the new_mean using the mean update formula.
  4. Calculate the new_M2 using the sum of squares update formula, leveraging old_mean and new_mean.
  5. Update mean to new_mean and M2 to new_M2 to prepare for the next data point.

Once all data points have been processed, the final mean is simply the last calculated mean. The variance can then be calculated from the final M2 and count.

• Population variance: ${\sigma }^2=M_k/k$
• Sample variance: $s^2=M_k/(k-1)$

As a data scientist, you'll almost always want the sample variance, as it provides an unbiased estimate of the true population variance from your sample of data. Note that this requires at least two data points (k > 1).

2. A Class-Based Implementation in Python

Given your background in software engineering, you'll recognize that this stateful, step-by-step process is perfectly suited for a class-based implementation. A class can encapsulate the state variables (count, mean, M2) and provide methods to update the state and retrieve the results.

The following resource provides an excellent, clear Python implementation of Welford's algorithm in a class named RunningStats.

Let's break down how the code in the RunningStats class maps directly to our formulas.

• self.n: This is our count, $k$.
• self.old_m: This is the mean before the update, ${\mu }_{k-1}$.
• self.new_m: This is the mean after the update, ${\mu }_k$.
• self.old_s: This is the sum of squared differences before the update, $M_{k-1}$.
• self.new_s: This is the sum of squared differences after the update, $M_k$.

Now, look at the push method in the else block:

# self.new_m = self.old_m + (x - self.old_m) / self.n
# self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m)

These two lines are a literal translation of the update formulas we derived. The implementation then updates old_m and old_s to prepare for the next call to push.

The variance() method correctly returns self.new_s / (self.n - 1), which is the sample variance, and handles the edge case where n < 2.

3. Refining the Implementation

The RunningStats class from the Stack Overflow post is perfectly functional, but it uses several variables (old_m, new_m, old_s, new_s) to manage the state within the push method. We can write this more concisely.

The Wikipedia article on this topic provides a slightly different, more compact implementation style.

Let's analyze the logic from the Wikipedia update function and see how it achieves the same result:

# pseudo-code adapted from Wikipedia
count += 1
delta = new_value - mean  # Here, `mean` is µ_{k-1}
mean += delta / count     # `mean` is now updated to µ_k
delta2 = new_value - mean # Here, `mean` is the new µ_k
M2 += delta * delta2

This sequence is clever. It calculates the difference (x_k - µ_{k-1}) and stores it as delta. It then updates the mean. Finally, it calculates the second difference (x_k - µ_k) as delta2 after the mean has been updated. This allows it to compute the product (x_k - µ_{k-1})(x_k - µ_k) without needing to store old_mean and new_mean separately.

We can refactor our push method to use this more efficient approach. Here is what a revised RunningStats class could look like:

import math

class Welford:
    """
    Implements Welford's online algorithm for computing mean and variance.
    """
    def __init__(self):
        self.count = 0
        self.mean = 0.0
        self.M2 = 0.0 # M2 is the sum of squares of differences from the current mean

    def update(self, new_value):
        self.count += 1
        delta = new_value - self.mean
        self.mean += delta / self.count
        delta2 = new_value - self.mean # Uses the *new* mean
        self.M2 += delta * delta2

    def get_mean(self):
        return self.mean

    def get_variance(self, sample=True):
        if self.count < 2:
            return 0.0 # Or float('nan')
        
        if sample:
            return self.M2 / (self.count - 1)
        else: # Population variance
            return self.M2 / self.count

    def get_stddev(self, sample=True):
        return math.sqrt(self.get_variance(sample))

This version is more memory-efficient and, for many, a more elegant expression of the algorithm. It is the version we will use as we move forward.

Conclusion

In this lesson, we successfully bridged the gap between theory and practice. We translated the mathematical recurrence relations for Welford's algorithm into a robust and reusable Python class, capable of computing mean and variance in a single pass over a data stream.

Key Takeaways:

• Welford's algorithm can be implemented by tracking three state variables: the count of items (count), the running mean (mean), and the sum of squared differences from the mean (M2).
• A class is an ideal structure for implementing online algorithms like Welford's, as it encapsulates the state and the logic for updating it.
• The update logic can be written concisely by calculating intermediate differences (delta and delta2) in a specific order, avoiding the need for separate "old" and "new" state variables within the update step.

Preview of the Next Lesson:

Now that we have a solid implementation of Welford's algorithm, our next step is to create the other algorithms we'll need for our final comparison. In the next lesson, you will implement the naive single-pass and the stable two-pass variance algorithms. This will complete our toolkit and set the stage for the final experiment where we will empirically demonstrate the superiority of Welford's algorithm on a carefully constructed dataset.