# Calculate variance from a stream of sample values

I'd like to calculate a standard deviation for a very large (but known) number of sample values, with the highest accuracy possible. The number of samples is larger than can be efficiently stored in memory.

The basic variance formula is:

$\sigma^2 = \frac{1}{N}\sum (x - \mu)^2$

... but this formulation depends on knowing the value of $\mu$ already.

$\mu$ can be calculated cumulatively -- that is, you can calculate the mean without storing every sample value. You just have to store their sum.

But to calculate the variance, is it necessary to store every sample value? Given a stream of samples, can I accumulate a calculation of the variance, without a need for memory of each sample? Put another way, is there a formulation of the variance which doesn't depend on foreknowledge of the exact value of $\mu$ before the whole sample set has been seen?

-
This very same issue was discussed on dsp.SE, and the answers were very similar, with numerically unstable methods proposed in one answer, and more stable methods described by others. – Dilip Sarwate Mar 4 '12 at 17:04

You can keep two running counters - one for $\sum_i x_i$ and another for $\sum_i x_i^2$. Since variance can be written as $$\sigma^2 = \frac{1}{N} \left[ \sum_i x_i^2 - \frac{(\sum_i x_i)^2}{N} \right]$$

you can compute the variance of the data that you have seen thus far with just these two counters. Note that the $N$ here is not the total length of all your samples but only the number of samples you have observed in the past.

-
Wait -- shouldn't the N be inside the the expected values, such that the right side is divided by $N^2$? – user6677 Feb 5 '11 at 22:05
e.g. $$\sigma^2 = \frac{(\sum_i x_i^2)}{N} - (\frac{\sum_i x_i}{N})^2$$ – user6677 Feb 5 '11 at 22:27
@user6677: You are indeed right. Thanks for the correction. – Dinesh Feb 5 '11 at 22:30
Note the sum of squares especially is at risk of overflow. – dfrankow Dec 29 '14 at 18:06
Instead of having counters for sum and sum of squares it would make more sense to have counters for average and average square. – Karolis Juodelė Feb 6 '15 at 5:50

I'm a little late to the party, but it appears that this method is pretty unstable, but that there is a method that allows for streaming computation of the variance without sacrificing numerical stability.

Cook describes a method from Knuth, the punchline of which is to initialize $m_1 = x_1$, and $v_1 = 0$, where $m_k$ is the mean of the first $k$ values. From there,

\begin{align*} m_k & = m_{k-1} + \frac{x_k - m_{k-1}}k \\ v_k & = v_{k-1} + (x_k - m_{k-1})(x_k - m_k) \end{align*}

The mean at this point is simply extracted as $m_k$, and the variance is $\sigma^2 = \frac{v_k}{k-1}$. It's easy to verify that it works for the mean, but I'm still working on grokking the variance.

-
+1, excellent! I didn't feel a simple upvote would be sufficient to express my appreciation of this answer, so there's an extra 50 rep bounty coming your way in 24 hours. – Ilmari Karonen Mar 4 '12 at 20:12
Isn't the final variance calculation should be Vk / k (not k-1)? – errr Oct 10 '13 at 15:49

tldr: I think there is a missing $(k-1)$ term in @Dan's variance calculation. I think it should be $$\sigma^2_k = \frac{1}{k} \left( (k-1) \sigma^2_{k-1} + (x_k-m_{k-1})(x_k-m_k)\right)$$

Validation: I tried what Dan's formula but am getting something different from Matlab's var command. For example, take 1000 Gaussian numbers (mean 0, var 100) and compute their mean and var:

clear; clc
n=1000;
x=round(randn(1,n)*10);

% calcs
mm(1)=x(1);
vm(1)=0;
m(1)=x(1);
v(1)=0;
ma(1)=x(1);
va(1)=0;
for k=2:n
mm(k) = mean(x(1:k)); %matlab mean
vm(k) = var(x(1:k),1); %matlab var

m(k) = m(k-1) + 1/k * (x(k)-m(k-1)) ; %Knuth mean
v(k) = (1/k) * (v(k-1) + (x(k)-m(k-1)) * (x(k)-m(k))) ; %Knuth var

ma(k) = (1-1/k) * ma(k-1) + (1/k) * x(k); %Finch mean
%va(k) = (1-1/k) * (va(k-1) + (1/k) * (x(k)-ma(k-1))^2 ); %Finch 143
va(k) = (1/k) * ((k-1)*va(k-1) + (x(k)-m(k-1)) * (x(k)-m(k)) ); %Finch 143 rearranged to look like Knuth
end

% plots
figure(1); clf
subplot 211
plot(1:n, x, 'k--')
hold on
plot(1:n,mm, 'b-')
plot(1:n,m, 'g-')
plot(1:n,ma, 'r--')
hold off
ylabel('means')
legend('x','Matlab','Knuth','Finch')

subplot 212
plot(1:n,vm, 'b-')
hold on
plot(1:n,v, 'g-')
plot(1:n,va, 'r--')
hold off
ylabel('vars')
xlabel('steps')
legend('Matlab','Knuth','Finch')


As you can see I selected Matlab's (octave's) VAR(X,1) command which divides by $N$ and not by $(N-1)$, which goes back to @errr's comment - it has to do with using biased vs unbiased estimator (I usually go to biased one because it makes more sense to me)

So it seems like Knuth variance is ok for a couple of steps and then it's way below because of the missing $(k-1)$ term.

So what I tried then is based on Tony Finch's explanation of recursive variance (eq (143)) with $\alpha=\frac{1}{k}$. Admittedly it's a hack as my alpha changes with every step, whereas his is assumed to be constant (?), but in the end this formula agrees with matlab's var(X,1).

-
If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. – Mark Fantini Jul 14 '14 at 22:56
I guess it was unclear - I found a mistake in previous answer. added a small summary to say that up front. Thank you! – alexey Jul 15 '14 at 0:01
@alexey I think your implementation of the Knuth recurrence does not match the one given in Dan Lecocq's answer. In his answer v(k) is not the variance; the variance can be computed from v(k). My tests using his recurrence are working fine. – rhashimoto Jan 30 at 22:27