# incrementally calculate gaussian distribution

Is there a way to incrementally calculate the parameters for a Gaussian distribution? The only parameters that are of interest are the mean and standard deviation, not the height (ie. normal distribution)

The problem is, I visit about 100 billion data points, of which I need to quickly calculate/update the distribution parameters, but I can not store all points in memory or visit all points for every estimation.

I'm looking for a way to calculate the distribution with only keeping a persistent state of size N(1) and the latest visited data point.

I can calculate the 'running mean' by this (assuming no numerical inaccuracies):

mean[0]=0;
mean[n+1] = mean[n] + (value - mean[n]) / (n+1)
#or alternatively: (more accurate in my experiments, but also slower)
mean[n+1] = (n * mean[n] + value) / (n+1)


but I don't know how to calculate the 'running deviation' (or variance for that matter)

I could build a histogram of values, but the problem is that the initial spread of the values is unknown and could span in the hundreds but could just as well span millions (depending on the specific dynamic system generating values). Furthermore, as the simulation goes on, the deviation typically increases monotonically, starting at 0.

-

Recall that the standard deviation $\sigma$ for some sample $X$ with results $x_{1},x_{2},\dots,x_{i}$ can be calculated by:

$$\sigma=\sqrt{\frac{\sum_{i}x_{i}^{2}}{n}-\left(\frac{\sum_{i}x_{i}}{n}\right)^{2}}=\sqrt{E(X^{2})-E^{2}(X)}$$

Therefore, if we add a new result $x_{\text{new}}$ to the sample, we can calculate the new mean $\bar{x}_{\text{new}}$ using:

$$\bar{x}_{\text{new}}=\frac{n\bar{x}+x_{\text{new}}}{n+1}$$

And then, rearranging $\sigma^{2}=E(X^{2})-E^{2}(X)=\frac{\sum_{i}x_{i}^{2}}{n}-\bar{x}^{2}$ gives $\sum_{i}x_{i}^{2}=n(\sigma^{2}+\bar{x}^{2})$, therefore, we find that:

$$\sigma^{2}_{\text{new}}=\frac{n(\sigma^{2}+\bar{x}^{2})+x_{\text{new}}^{2}}{n+1}-\bar{x}_{\text{new}}^{2}$$

Or, in pseudocode:

mean[n+1] = ((n * mean[n]) + value)/(n+1)
var[n+1] = (n * (var[n] + mean[n]^2) + value^2)/(n+1) - mean[n+1]^2

-
Thanks. Though I don't recall that formula. I recall $$\sigma=\sqrt{\frac{\sum_{i}\left(x_{i}-\bar{x}\right)^2}{n}}$$. I also fixed your pseudocode which wasn't the same as the formula. –  Mark Jeronimus Feb 17 '13 at 14:37
@Zom-B The two formulae are equivalent. And thank you for fixing my pseudocode, I missed that!! –  Shaktal Feb 17 '13 at 14:42
I'm now at this stage. I never thought it would turn out this small. For each sample: $$Sum_{n+1}=Sum_{n}+value$$$$SqSum_{n+1}=SqSum_{n}+value^2$$ and getting the parameters: $$mean=\frac{Sum}{n}$$$$stdev=\sqrt{\frac{SqSum_{n}}{n} - (\frac{Sum_{n}}{n})^2}$$ –  Mark Jeronimus Feb 17 '13 at 15:51

Here is an updating formula for the variance $\sigma_n^2$ (i.e., the square of the SD) in terms of the mean $\mu_n$:

$$(n+1)\sigma_{n+1}^2 = n \sigma_n^2 + n \left ( \frac{\mu_n}{n+1} + x_{n+1} \right )^2 + (x_{n+1}-\mu_{n+1})^2$$

This is derived by considering

$$n \sigma_n^2 = (x_1 - \mu_n)^2 + \ldots + (x_n - \mu_n)^2$$

\begin{align}(n+1) \sigma_{n+1}^2 &= (x_1 - \mu_n)^2 + \ldots + (x_n - \mu_n)^2+ (x_{n+1} - \mu_{n+1})^2\\ &= \left ( x_1 - \frac{n \mu_n + x_{n+1}}{n+1} \right )^2 + \ldots + (x_{n+1} - \mu_{n+1})^2\\ &= (x_1-\mu_n)^2 + 2 (x_1-\mu_n) \left (\frac{\mu_n}{n+1} + x_{n+1} \right ) + \left (\frac{\mu_n}{n+1} + x_{n+1} \right )^2 + \dots+ (x_{n+1} - \mu_{n+1})^2\\ &= n \sigma_n^2 + n \left(\frac{\mu_n}{n+1} + x_{n+1} \right )^2 + (x_{n+1} - \mu_{n+1})^2 \\ \end{align}

because $\sum_{k=1}^n (x_k-\mu_n) = 0$.

In code:

nvar[n+1] = nvar[n] + n*((mean[n]/(n+1) + value)^2 + (value-mean[n+1])^2;
sd[n+1]=sqrt(nvar[n+1]/(n+1));

-
Your pseudocode and formula are not equivalent and I can get neither to work. As Shaktal's formula is simpler I'm not going to bother deriving it again the way you did. –  Mark Jeronimus Feb 17 '13 at 14:51
Sorry for the inconvenience, I left off the square in the second term. Otherwise it is fine. But feel free to ignore as you need. Interesting problem, though. –  Ron Gordon Feb 17 '13 at 14:55