01Brief10 / 14

Dividing by n, or by n minus one

Write spread(values, sample=True) returning the mean, variance and standard deviation, dividing by n − 1 for a sample and by n for a whole population.

Function	spread()
Track	Data without the library
Difficulty
Time	~22 min
Concepts	Bessel’s correction sample versus population two passes over the data
Cases	11

02The contract

Return a dict with mean, variance and deviation.
The variance is the mean of the squared distances from the mean.
For a sample — the default — divide by n − 1, not n.
For a whole population, divide by n.
The deviation is the square root of the variance.
A single value has a population variance of 0 and no sample variance at all: raise ValueError.
No values at all: raise ValueError.
Leave the caller’s list exactly as you found it.

03Why this one

Variance is the first formula where the textbook and the library disagree with each other and both are right. numpy.var divides by n; statistics.variance and pandas.Series.var divide by n − 1; R divides by n − 1. The difference is not pedantry — it is the difference between describing the numbers you have and estimating the population they came from, and on small samples it changes the answer materially. Every data scientist meets this, usually by noticing two tools disagree.

—The concept behind it

Lists on pythonmastery.io →

That is the reading half. Read it there if the idea is new, then come back here — explanation is what makes something make sense, and practice is what makes it stay.

—Worth doing first

The average, and the other average

Not a gate — start this one whenever you like. These cover ideas it leans on, and a tick means you have already done them.

—What this leads to

How unusual is this number

These build on what you are about to write here.

04Study first: why the sample divides by one fewer

The same four numbers, both ways. The gap is a quarter of the answer, and it grows as the sample shrinks.

values = [2, 4, 4, 6]
mean = sum(values) / len(values)          # -> 4.0
squares = [(v - mean) ** 2 for v in values]   # -> [4.0, 0.0, 0.0, 4.0]

sum(squares) / len(values)        # -> 2.0     population
sum(squares) / (len(values) - 1)  # -> 2.667   sample

# The mean was computed from the same four numbers, so the distances
# from it are systematically a little too small. Dividing by n - 1
# corrects for that.
#
# With one value: sum(squares) is 0 and n - 1 is 0. Nothing to divide.

The mean is an estimate, and the data already fitted it
The sample mean is the value that makes the sum of squared distances as small as it can possibly be — for this sample. So the spread you measure around it is biased downwards as an estimate of the real population spread. Dividing by n − 1 rather than n scales it back up by exactly the right amount. That is Bessel’s correction, and 'we used up one degree of freedom estimating the mean' is the short version.
Which one you want depends on what the numbers are
If your list *is* the whole thing — every pixel in this image, every employee in this company — divide by n, because there is no population beyond it to estimate. If it is a sample standing in for something larger, divide by n − 1. numpy defaults to n and statistics defaults to n − 1, which is why two tools can disagree about the same data and neither is broken.
One value has no sample variance
n − 1 is zero, and the honest answer is not 0 — it is that a single observation tells you nothing about spread. Returning 0 would claim a population is perfectly uniform on the evidence of one member. The population variance of one value genuinely is 0, because that one value is the whole population, which is a nice demonstration that the two formulas answer different questions.

your_code.py

05One way to write it

import math


def spread(values, sample=True):
    count = len(values)
    divisor = count - 1 if sample else count

    # Covers both refusals at once: an empty list either way, and a
    # single value when it is meant to be a sample.
    if divisor <= 0:
        raise ValueError(
            f"not enough values to measure spread: {count} "
            f"{'sample' if sample else 'population'} value(s)"
        )

    mean = sum(values) / count

    # Two passes, deliberately: the one-pass formula is the same algebra
    # and much worse arithmetic.
    total = sum((value - mean) ** 2 for value in values)
    variance = total / divisor

    return {
        "mean": mean,
        "variance": variance,
        "deviation": math.sqrt(variance),
    }

One guard covers both refusals

`divisor <= 0` is true for an empty list in either mode and for a single value in sample mode — which is exactly the set of inputs with no answer. Computing the divisor before validating is what makes that possible, and it means the two contract lines about refusing are one line of code that cannot drift apart.

The mean always divides by count, never by the divisor

Bessel’s correction applies to the variance and not to the mean — the mean of four numbers is their total over four whatever you intend to do next. Using the divisor for both is an easy slip once the variable is sitting there, and it produces a mean that is wrong by a factor of n/(n−1).

Two passes beats the clever one-pass formula

E[X²] − E[X]² is the same quantity and catastrophically unstable: on values like 100000001 and 100000002 it subtracts two nearly equal large numbers and the answer can come out negative, which a variance cannot be. Welford’s algorithm is the genuinely good single-pass method if you need streaming. Never the naive one.

06When not to write this at all

statistics.variance and statistics.pvariance are the two functions here, already written, and statistics.stdev and pstdev give the deviations — the names encoding which divisor each uses, with the p for population. numpy.var and pandas both take a ddof argument instead, defaulting to 0 and 1 respectively, which is the single most common cause of two pipelines disagreeing about the same column. Whichever you use, write down which convention you meant. And if the data arrives as a stream too large to hold, look up Welford rather than trying to make the two-pass version fit.

07Take it further

Implement Welford’s single-pass algorithm and check it agrees with this one to many decimal places on a few thousand values.
Show the naive one-pass formula failing: feed both versions [100000001, 100000002, 100000003] and look at what comes out.
Add covariance between two lists, and notice which divisor it needs and why it is the same argument.