Streaming statistics

Situations where streaming statistics are useful:

For the remainder, let’s consider a set of observations yiy_i, weights wiw_i, such that

 0 

A naive approach to calculating a weighted streaming mean, μ^\widehat{\mu} and unbiased streaming variance, , would be to calculate:

 1 

where

 2 

This calculation however does not hold for large nn values.

An alternative calculation was suggested by West1 cite:&west1979updating West, D. (1979). Updating mean and variance estimates: an improved method. Communications of the ACM, 22(9), 532–535. , where we calculate:

 3 

Let’s look at an example in Python.

from plotutils import *
import numpy as np

mu = 10.0
sigma = 20.0
N = 100_000

Y = np.random.normal(loc=mu, scale=sigma, size=N)
import matplotlib.pyplot as plt
import matplotlib

matplotlib.rc('figure', figsize=(12, 6))
plt.hist(Y, bins=50,color="lightpink")
plt.show()

As expected, the mean and variance when calculated in “batch” mode should be equivalent to the original values μ\mu and σ\sigma.

print(f"mean: {np.mean(Y)}")
print(f"std: {np.std(Y)}")
mean: 10.040100563607174
std: 19.946298541957912
class StreamingStatistcs:
    def __init__(self):
        self.sum = 0.0
        self.mean = 0.0
        self.t = 0
        self.n = 0
        self.var = None

    def calculate(self, y, w):
        q = y - self.mean
        tmp_sum = self.sum + w
        r = q*w / tmp_sum

        self.mean += r
        self.t += q*r*self.sum
        self.sum = tmp_sum
        self.n += 1
        if (self.sum == 0.0 or self.n < 2):
            self.var = 0
        else:
            self.var = (self.t*self.n)/(self.sum*(self.n-1))
        return (self.mean, self.var)

means = []
vars = []
st = StreamingStatistcs()
for y in Y:
    m, v = st.calculate(y, 1.0)
    means.append(m)
    vars.append(np.sqrt(v))
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.suptitle('Streaming statistics for $N=10^5$ observations')
ax1.plot(means)
ax1.hlines(y=10, xmin=0, xmax=N, colors="red")
ax1.set_title("Streaming mean")
ax2.plot(vars)
ax2.hlines(y=20, xmin=0, xmax=N, colors="red")
ax2.set_title("Streaming variance")
plt.show()