We calculate descriptive statistics—count, mean, unbiased variance—both overall and in a sliding window, of sequences of numerical data. The statistics calculator is written as a foldable accumulator function, decoupling calculation from iteration. Exactly the same stateless accumulator function works over data distributed over space in memory or in tables, or over time in asynchronous or lazy streams. Foldable accumulator functions can be tested independently of data source, say over ground-truth test data, then later deployed in harsh asynchronous environments with high assurance, often without even being recompiled.
The solution at the end of this notebook has been implemented in C, Python, and Clojure and deployed in mission-critical web applications and embedded applications. It runs in
O(w)
space, where
w
is the width of the window, and in
O(1)
time. Contrast to more obvious solutions that run in
O(w)
time.

Prelude: Running Count

For ground truth, we want a fixed but random sequence of data. The following is a sample of length 10 from the standard normal distribution—mean zero, standard deviation one:
In[]:=
zs={-0.178654,0.828305,0.0592247,-0.0121089,-1.48014,-0.315044,-0.324796,-0.676357,0.16301,-0.858164};
Fold has three arguments: an accumulator function, an initial state, and a sequence of data. All variations of Fold—over asynchronous streams, over on-demand (lazy) streams, over tables, lists, arrays—have exactly the same signature as Fold.
In[]:=
ClearAll[cume];​​cume[n_,z_]:=n+1;​​Fold[cume,0,zs]
Out[]=
10
The first argument of the accumulator function, cume, is a circulating state, kept externally to the function. Above, the circulating state is
n
. The caller must feed in two arguments to cume: the prior value
n
of the circulating state and a new observation
z
. The caller receives the posterior value of the state,
n+1
.
To see all intermediate results, use FoldList:
In[]:=
FoldList[cume,0,zs]
Out[]=
{0,1,2,3,4,5,6,7,8,9,10}

LEFT FOLD

Wolfram’s documentation for Fold includes an example that explains this operator symbolically. The following should be self-explanatory:
In[]:=
Fold[f,z,{a,b,c,d}]
Out[]=
f[f[f[f[z,a],b],c],d]
Fold first computes
f[z,a]
, then feeds the result—the posterior value of the circulating state—back into
f
, computing
f[f[z,a],b]
, and so on. More concisely, Fold iteratively applies the binary function
f
to the sequence
{a,b,c,d}
. That’s why we say Fold decouples calculation from iteration. Fold places no restriction on its accumulator-function argument
f
other than that it accept an argument of the type
T
z
of
z
in its first position, an argument of the type
T
a
of any of
a,b,c,d
, in its second position, and produces a result of type
T
a
.

RIGHT FOLD

It is important to know that there is a dual: the right fold. Left and Right folds can behave differently on lazy lists, so be on the lookout! We use only left fold in this document.
In[]:=
Fold[{x,y}f[y,x],z,{a,b,c,d}]
Out[]=
f[d,f[c,f[b,f[a,z]]]]

Running Mean

Write
x
for the mean-so-far, and write a new accumulator function that computes the running mean. Its circulating state includes the running mean
x
, the running count
n
, and the running sum. This new accumulator function computes the new values of all three state variables in the obvious way. Note that all three state variables are necessary to track the running mean calculated this way. We relieve that situation shortly below.
In[]:=
ClearAll[cume];​​cume[{x_,n_,sum_},z_]:=​​
sum+z
n+1
,n+1,sum+z;​​Fold[cume,{0,0,0},zs]
Out[]=
{-0.279472,10,-2.79472}
Check against Mathematica built-in.
In[]:=
Mean[zs]
Out[]=
-0.279472
Notice how the mean stabilizes as data accumulate, here over a longer and random sample:
In[]:=
ListLinePlot[#〚1〛&/@FoldList[cume,{0,0,0},RandomVariate[NormalDistribution[],10000]]]
Out[]=

Preferable Form

Here is a preferable form with one fewer quantity in the circulating state. This form is so powerful that we state it as a mantra or slogan:
compute the new value as sum of the old value and a correction that depends only on old values and the single observation
z
.
We will strive for this form in everything that follows. It is the essential feature of all recurrent filter algorithms.
The right-hand side of the recurrence
z
←
z
+K×(z-
z
)
,
z
here standing for the running mean, depends only on old values and the new observation
z
. The recurrence is not an equation, but it’s shorthand the following equation:
z
n+1
=
z
n
+K×(
z
n+1
-
z
n
)=x+K×(
z
n+1
-x).
(
1
)
In[]:=
ClearAll[cume];​​cume[{x_,n_},z_]:=​​WithK=
1
n+1
,​​{x+K(z-x),n+1};​​Fold[cume,{0,0},zs]
Out[]=
{-0.279472,10}
We prefer this form because (1) it separates old quantities on the right (except for the new datum,
z
) from new quantities on the left, and (2) it is easy to memorize because it is an affine update: the old estimate plus a correction written as a gain
K
times the residual
z-
z
. The residual is the difference between the new observation and the old mean-so-far.
To prove this formula
z
+K×(z-
z
)
correct, we need only show that it equals
sum+z
n+1
, which is obviously the new mean. Write
z
, the old mean, as
z
(n+1)
n+1
=
n
z
+
z
n+1
and write
z
+K×(z-
z
)=
n
z
+
z
n+1
+
z-
z
n+1
=
sum+z
n+1
, noting that the sum is
n
z
. I think this is beautiful. The same reasoning will help us derive many more sophisticated formulas, up-to-and-including the Kalman update.
For another twist on the proof above, we can let Mathematica do the work. The Kalman-like form
x+K(z-x)
is equal to the obvious form
(nx+z)/(n+1)
when
n+1≠0
.
In[]:=
Solvex+K(z-x)
nx+z
n+1
,K
Out[]=
K
1
1+n


Running Variance

Definitions

Variance is the sum of squared residuals against the running mean, divided by the count less one. The “count less one” is Bessel's correction.
In[]:=
Variance[zs]
Out[]=
0.395183
Variance is also the square of standard deviation.
In[]:=
2
StandardDeviation[zs]
Variance[zs]
Out[]=
True
The sum of squared residuals,
Σ
n
def
=
n
∑
i=1
2
(
z
i
-
z
n
)
(
2
)
is the variance times the length-less-one:
In[]:=
Variance[zs]*(Length[zs]-1)
Out[]=
3.55665
The sum of squared residuals will be a cornerstone quantity for all future calculations. Here is another look at it.
In[]:=
With[{μ=Mean[zs]},​​Plus@@(
2
(#-μ)
&/@zs)]
Out[]=
3.55665

Requirement

We require an incremental, running variance that does not refer to future observations—even indirectly, through the eventual mean
z
n
—and does not store the entire past, only a summary of constant size.
The definition explicitly refers to the past history. It seems, at first glance, impossible to compute the variance without keeping the entire history. But it possible, and we’ll get there in baby steps.

Brute-Force Variance (doesn’t meet requirements)

The following accumulates variance directly from the definition, useful for numerical checks. However, it assumes we know length and mean of the entire sequence zs ahead of time, so it does not meet the requirement of being incremental. It also needs computer memory for the entire sequence zs, which is of variable, unknown, or infinite length, so does not meet the requirement of constant memory.
◼
  • Re-use the variable zs, now for a shorter sequence of ground-truth data, convenient for checking arithmetic in one’s head.
    • In[]:=
    ClearAll[zs];zs={55,89,144};
    In[]:=
    Variance[zs]//N
    Out[]=
    2017.
    Hand-calculated mean:
    In[]:=
    (55+89+144)/3.
    Out[]=
    96.
    Hand-calculated variance = SSR (sum of squared residuals)
    /(n-1)
    :
    In[]:=
    ((55-96.)^2+(89-96.)^2+(144-96.)^2)/2.
    Out[]=
    2017.
    Now, we’ll write a truly nasty accumulator function to produce the final variance. It’s nasty because it refers to the entire sequence through the variable zs.
    In[]:=
    ClearAll[cume];​​cume[var_,z_]:=var+
    2
    (z-Mean[zs])
    Length[zs]-1
    ;​​FoldList[cume,0.,zs]
    Out[]=
    {0.,840.5,865.,2017.}
    It’s even nastier because the intermediate variances are not valid. They do not refer to the “mean-so-far”—the incremental running mean, but to the final mean.

    School Variance (meets requirements)

    The following is much better, exploiting the “school formula,” proved by expanding the square. Letting
    z
    n
    =x
    , incremental running mean, the following yields the sum of squared residuals, ssq:
    n
    ∑
    i=1
    2
    (
    z
    i
    -
    z
    n
    )
    =
    n
    ∑
    i=1
    2
    z
    i
    -n
    2
    z
    n
    (
    3
    )
    The variance computed this way is incremental and runs in constant memory, though it is not yet a recurrence, i.e., old value plus a correction.

    Recurrent Variance (meets requirements)

    Symbolic Proof

    Here’s the old ssq from the School formula:
    Here’s the new ssq via the same formula:
    Transcribe Equation 6.
    Transcribe the “obvious equalities” above.
    Note that the rules don’t do much to the definition of s[n].
    The rules instantly yield the correction terms in Recurrence 4 and Equation 5.

    Proof by Hand

    Remember the “obvious equalities.”
    Expand the second term:
    Cancel terms and complete the square:

    Equivalence to Welford’s

    Welford’s famous recurrence is
    We won’t write out the algebraic manipulations, but it’s good practice.

    Folding It

    Recurrence 4 lets us get rid of one level of With, one level of dependency:

    Windowed Statistics

    Foldable

    runningStatsFoldable is a a foldable accumulator function that effects windowed statistics. Its circulating state carries extra information for pedagogy and convenience. It’s possible to make the interface much tighter. The derivation of its formulae are written in comments. The specific example data in this section are taken from other deployments of this code in C and Python. In the next section, we add some checks of the results against Mathematica’s built-ins.

    Random Data, Checked

    The solution does not work when the window size is zero, but is robust even in the nonsensical case of the window’s being longer than the data.

    CITE THIS NOTEBOOK

    Windowed incremental online statistics: extracting models from data, one observation at a time​
    by Brian Beckman​
    Wolfram Community, STAFF PICKS, May 12, 2024
    ​https://community.wolfram.com/groups/-/m/t/3174136