Standard deviation for online scenarios

Question

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.

Answer 1

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

    def add_sample(self, x):
        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