# Standard deviation for online scenarios

Let's suppose you have to create an online algorithm $$f(x_i)$$ that will return

• the average ($$\mu$$); And
• the standard deviation ($$\sigma$$);

Of numbers $$x_0, \dots, x_i$$.

A straightforward approach would be to simply store the argument $$x_i$$ into an array ($$a$$), update a global counter ($$c$$) to recalculate $$\mu$$, and then recalculate $$\sigma$$. As follows:

$$f(x_i):$$

1. $$a \gets a \cup \{x_i\}$$;
2. $$c \gets c + x_i$$;
3. $$\mu = \frac{c}{i}$$;
4. $$\sigma = \sqrt{\frac{1}{i} \sum_{k = 0}^i (a_k - \mu)^2}$$;
5. return $$\mu$$, and $$\sigma$$;

The time complexity of the above algorithm in an iteration $$i$$ is $$O(i)$$, and the "global" space complexity is $$O(n)$$.

I was wondering if there are strategies for reducing the time and "global" space complexities. Precisely, if it is possible to simply disconsider the array $$a$$, and to avoid iterating over all the previous numbers for recalculating the STD.

There are "streaming" online versions for computing at least the first four discrete statistical moments (mean, variance, skewness, kurtosis)

The mean is relatively easy to see how it operates:

1. Store the current mean $$M$$ and sample count $$N$$.
2. To add a sample $$v$$, $$\mathrm{M} \gets (\mathrm{M} N + v)/(N+1)$$, $$N \gets N + 1$$

For more details on the higher moments, see this link.

Here is a Python implementation which computes these first 4 moments in an online manner (with sample variants):

class stream_stats:
def __init__(self):
self.n = 0
self.mean = 0
self.M2 = 0
self.M3 = 0
self.M4 = 0

n1 = self.n
self.n = self.n + 1
delta = x - self.mean
delta_n = delta / self.n
delta_n2 = delta_n**2
term1 = delta * delta_n * n1
self.mean = self.mean + delta_n
self.M4 = (
self.M4
+ term1 * delta_n2 * (self.n**2 - 3 * self.n + 3)
+ 6 * delta_n2 * self.M2
- 4 * delta_n * self.M3
)
self.M3 = self.M3 + term1 * delta_n * (self.n - 2) - 3 * delta_n * self.M2
self.M2 = self.M2 + term1

def variance(self, DOF=1):
return self.M2 / (self.n - DOF)

def stddev(self, DOF=1):
from math import sqrt

return sqrt(self.variance(DOF))

def skewness(self, sample=True):
# Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
if sample:
return (
self.n
* self.M3
/ ((self.n - 1) * (self.n - 2) * (self.M2 / (self.n - 1)) ** 1.5)
)
else:
return self.M3 / (self.n * (self.M2 / self.n) ** 1.5)

def kurtosis(self, sample=True):
# Caution: If all the inputs are the same, M2 will be 0, resulting in a division by 0.
if sample:
return (
(self.n + 1)
* self.n
* (self.n - 1)
* self.M4
/ ((self.n - 2) * (self.n - 3) * self.M2**2)
)
else:
return (self.n * self.M4) / (self.M2**2) - 3

def excess_kurtosis(self, sample=True):
if sample:
return self.kurtosis(sample) - 3 * (self.n - 1) ** 2 / (
(self.n - 2) * (self.n - 3)
)
else:
return self.kurtosis(sample) - 3

• If I may add one more question. Have you played with such algorithm with C/C++? I've playing with a C++ implementation (johndcook.com/blog/skewness_kurtosis) and I am obtaining, almost, right results for the sample [7,8,0,4,3,6,4,1,0.5]. It seems it is related to the numerical stability issue. Apr 19 at 14:06
• Oh, I found a bug in the non-sampling implementation of how I computed skewness. Fixed it. As far as numerical accuracy, when I test your sequence it's within a few machine epsilon of what scipy.stats reports. This implementation is just a port of what I use in another C++ code I use. Apr 19 at 18:13