Comparing Single-Pass and Two-Pass Variance Algorithms

1. Implementing the "Naive" Single-Pass Algorithm

First, let's implement the algorithm that is most directly derived from the "textbook" formula for variance. As we discussed in the first module, this approach is fast but notoriously prone to numerical errors due to catastrophic cancellation.

The formula, often called the "computational formula" for sample variance, is:
$s^2=\frac{{\sum }_{i=1}^n{x}_i^2-\frac{({\sum }_{i=1}^n{x}_i{)}^2}{n}}{n-1}$

To implement this, we need to iterate through the data just once, keeping track of three values:

1. The number of data points, $n$.
2. The sum of the data points, $\sum x_i$.
3. The sum of the squares of the data points, $\sum x_i^2$.

After the loop, we plug these three aggregates into the formula. The Wikipedia article on this topic provides excellent pseudocode for this exact algorithm.

Following that structure, here is a straightforward Python implementation. This function takes a list of numbers and returns the sample variance.

def naive_variance(data, sample=True):
    """
    Calculates variance using the naive, single-pass algorithm.
    This method is prone to numerical instability.
    
    Args:
        data (list or iterable): A list of numerical data.
        sample (bool): If True, calculates sample variance (n-1 denominator). 
                       Otherwise, calculates population variance (n denominator).

    Returns:
        float: The calculated variance.
    """
    n = 0
    sum_val = 0.0
    sum_sq = 0.0

    for x in data:
        n += 1
        sum_val += x
        sum_sq += x**2

    if n < 2:
        return 0.0 # Or float('nan')

    # Using the computational formula for variance
    if sample:
        variance = (sum_sq - (sum_val**2) / n) / (n - 1)
    else:
        variance = (sum_sq - (sum_val**2) / n) / n
        
    return variance

This implementation is clean and efficient in terms of passes, but as we know, its reliance on the subtraction of two potentially very large and similar numbers (sum_sq and (sum_val**2) / n) is its fatal flaw.

2. Implementing the Two-Pass Algorithm

Next, we'll implement the two-pass algorithm. This method is a significant improvement in numerical stability over the naive approach because it avoids the large-number subtraction problem. It adheres more closely to the definitional formula for variance.

The process involves two distinct passes over the data:

1. First Pass: Calculate the mean, $\mu =\frac{\sum x_i}{n}$.
2. Second Pass: Calculate the sum of squared differences from that mean, $\sum (x_i-\mu {)}^2$.

The sample variance is then this sum divided by $n-1$.

The primary drawback, of course, is the need to iterate over the data twice. This is inefficient if the dataset is massive and stored on slow media, and it's impossible for data streams that you can't store. However, for moderately sized datasets in memory, it serves as an excellent, stable baseline for comparison.

The same Wikipedia resource also provides a clear implementation of this algorithm.

Based on that logic, here is the Python function for the two-pass algorithm:

def two_pass_variance(data, sample=True):
    """
    Calculates variance using the stable, two-pass algorithm.

    Args:
        data (list or iterable): A list of numerical data.
        sample (bool): If True, calculates sample variance (n-1 denominator). 
                       Otherwise, calculates population variance (n denominator).

    Returns:
        float: The calculated variance.
    """
    n = len(data)
    if n < 2:
        return 0.0 # Or float('nan')

    # Pass 1: Calculate the mean
    mean = sum(data) / n

    # Pass 2: Calculate the sum of squared differences
    sum_sq_diff = sum((x - mean)**2 for x in data)

    # Final calculation
    if sample:
        variance = sum_sq_diff / (n - 1)
    else:
        variance = sum_sq_diff / n
        
    return variance

This function is more robust than our naive implementation and will be our "ground truth" or stable baseline when we later compare results.

Conclusion

In this lesson, we have successfully implemented the two key algorithms for comparison against Welford's method. You now have three distinct functions ready for our empirical analysis.

Key Takeaways:

• The naive single-pass algorithm is implemented by accumulating the count, sum, and sum of squares, then applying the computational formula. It's fast but numerically unstable.
• The two-pass algorithm is implemented by first computing the mean, then iterating again to compute the sum of squared differences. It's numerically stable but requires two passes over the data.
• You now have a complete toolkit: the Welford class from the previous lesson, and the naive_variance and two_pass_variance functions from this one.

Preview of the Next Lesson:

With our algorithmic toolkit complete, we can now move to the most interesting part: demonstrating the failure of the naive algorithm. In the next lesson, we will construct a synthetic dataset specifically designed to cause catastrophic cancellation. This will set the stage for the final lesson, where we run all three algorithms on this dataset and compare their accuracy.