Online estimation of variance with limited memory

Question

I am creating a component that aims to calculate the average and variance of a metric associated with events happening during time but with a limited internal memory.

Imagine that the events are visitors entering in a shop and the metric is their age.

During time, my component receives events with the age of each visitor. I don't want my component to memorize the history of each ages. Ideally, I would like a light component storing only: the average A, the variance V and the number of events N.

After each event with age E, I want to update those three values :

N<=N+1
A<=(A*N+E)/(N+1)
V<=???

What for V? I am thinking of something like :

V<=(V*N+(E-A)^2)/(N+1)

I know it is not exact as my previous V is using the old A which is no more the average.

Q1 - Is there an exact formula?
Q2 - If not, is my proposal a good estimate? Is it biased? Will it converge correctly when N increases?
Q3 - Is there a better formula?

Answer 1

Nice and simple algorithm for computing variance in online manner was described by Welford (1962). Below you can see C++/Rcpp implementation of it that works offline, but can be easily adapted to online scenario:

List welford_cpp(NumericVector x) {
  
  int n = x.length();
  double delta;
  double msq = 0;
  double mean = x[0];
  
  if (n > 1) {
    for (int i = 1; i < n; i++) { 
      delta = x[i] - mean;
      mean += delta / (i+1);
      msq += delta * (x[i] - mean);
    }
    return Rcpp::List::create(Rcpp::Named("mean") = mean,
                              Rcpp::Named("variance") = msq / (n-1));
  }
  
  return Rcpp::List::create(Rcpp::Named("mean") = mean,
                            Rcpp::Named("variance") = NAN);
}

As you can see, it needs to store only four variables: n, delta, msq and mean and computes mean and variance simultaneously as you wanted.

---

Welford, B. P. (1962). Note on a method for calculating corrected sums of squares and products. Technometrics 4(3): 419-420.

Answer 2

The variance can be expressed as proportional to the squared difference between every value and the mean value, or (as many threads here in stats.SE documented, like this answer I wrote to another question) it can alternatively be expressed as proportional to the squared pairwise difference between every sample.

So we know:

$$\text {Var}(x) = \frac{1}{n} \cdot \sum_{i}(X_i-\overline{X})^2 = \frac{1}{2n^2} \cdot \sum_{i,j}(X_i-X_j)^2$$

Let's say you add another sample, indexed as the last index, $k$. Your previous variance would be:

$$\text {Var}_{old}(x) = \frac{1}{2(n-1)^2} \cdot \sum_{i<k,j<k}(X_i-X_j)^2$$

Your new variance is

$$\text {Var}_{new}(x) = \frac{1}{2n^2} \cdot \sum_{i,j}(X_i-X_j)^2 = \frac{1}{2n^2} \cdot \left( \sum_{i<k,j<k}(X_i-X_j)^2 + \sum_{j<k}(X_k-X_j)^2 + \sum_{i<k}(X_i-X_k)^2 \right)$$

But

$$\sum_{j<k}(X_k-X_j)^2 = \sum_{i<k}(X_i-X_k)^2\\ \sum_{i<k,j<k}(X_i-X_j)^2 = {2(n-1)^2} \cdot \text {Var}_{old}(x) $$

So

$$\text {Var}_{new}(x) = \left(\frac{n-1}{n}\right)^2 \text{Var}_{old}(x)+ \frac{1}{n^2} \sum_{j<k}(X_k-X_j)^2 $$

As @MarkL.Stone said in the comments, this still isn't efficient because we must keep every $X_i$. So, let's expand the formula to arrive at something more tractable.

$$\frac{1}{n^2}\sum_{j<k}(X_k-X_j)^2=\\ \frac{1}{n^2}\sum_{j<k}(X_k^2 - 2 \cdot X_j \cdot X_k + X_j^2)=\\ \frac{1}{n^2}\left(\sum_{j<k}X_k^2 - 2 \cdot X_k \cdot \sum_{j<k} X_j + \sum_{j<k} X_j^2\right)=\\ \frac{1}{n^2}\left(k\cdot X_k^2 - 2 \cdot X_k \cdot (k-1) \cdot \overline{X_{old}} + (k-1) \cdot \overline{X^2_{old}}\right)\\ $$ Because $$\sum_{j<k} X_j = (k-1) \cdot \overline{X_{old}}\\ \sum_{j<k} X_j^2 = (k-1) \cdot \overline{X^2_{old}}$$

The final form is then

$$\text {Var}_{new}(x) = \left(\frac{n-1}{n}\right)^2 \text{Var}_{old}(x)+ \frac{1}{n^2}\left(k\cdot X_k^2 - 2 \cdot X_k \cdot (k-1) \cdot \overline{X_{old}} + (k-1) \cdot \overline{X^2_{old}}\right) $$

You can use this formula to update the variance effectively memory-wise. You can also complement it to use batches instead of single point updates.

Basically you need to store the average, the average of the squared samples, and the variance every iteration, and use it to update the variance formula.

---

Further

$$\overline{X^2_{old}} = \text{Var}_{old}(x) + (\overline{X_{old}})^2\\ \therefore\text {Var}_{new}(x) = \left(\frac{n-1}{n}\right)^2 \text{Var}_{old}(x)+ \frac{1}{n^2}\left(k\cdot X_k^2 - 2 \cdot X_k \cdot (k-1) \cdot \overline{X_{old}} + (k-1) \cdot \left(\text{Var}_{old}(x) + (\overline{X_{old}})^2\right)\right) $$

Which brings the number of quantities that need to be stored down to 2.

Answer 3

OK Andy W gave the answer. By conserving the $E^2$ average in the same way as the E average, you can use $V = exp(E^2)-exp(E)^2$.