Welford’s Algorithm

Computation of Sample Mean and Variance

In HFT, the first and second-order statistics often need to be determined. Necessarily, the algorithm that performs the computation must be O(1) and numerically stable. The sample mean and variance of a signal $p$ are signals given by $$\begin{array}{ccc}{\mu }_n& =& \frac{1}{n}\sum _{i=\mathrm{1..}n}{p}_i\\ {\sigma }_n^2& =& \frac{1}{n}\sum _{i=\mathrm{1..}n}{\left(p_i–{\mu }_n\right)}^2\end{array}$$Let $S_n=n{\mu }_n$ and $M_n=n{\sigma }_n^2$. Then $$\begin{array}{ccc}{S}_n–S_{n–1}& =& p_n\\ {\mu }_n& =& \frac{n–1}{n}{\mu }_{n–1}+\frac{1}{n}{p}_n\\ & =& {\mu }_{n–1}+\frac{1}{n}\left(p_n–{\mu }_{n–1}\right)\end{array}$$Also, $$\begin{array}{ccc}{M}_n–M_{n–1}& =& \sum _{i=\mathrm{1..}n}(p_i–{\mu }_n{)}^2–\sum _{i=\mathrm{1..}n–1}(p_i–{\mu }_{n–1}{)}^2\\ & =& \sum _{i=\mathrm{1..}n–1}\left[(p_i–{\mu }_n{)}^2–(p_i–{\mu }_{n–1}{)}^2\right]+(p_n–{\mu }_n{)}^2\\ & =& \sum _{i=\mathrm{1..}n–1}\left(2p_i–{\mu }_n–{\mu }_{n–1}\right)\left({\mu }_{n–1}–{\mu }_n\right)+(p_n–{\mu }_n{)}^2\\ & =& \left(n–1\right){\left({\mu }_{n–1}–{\mu }_n\right)}^2+(p_n–{\mu }_n{)}^2\\ & =& \frac{n–1}{n}\left(p_n–{\mu }_{n–1}\right)\left({\mu }_n–{\mu }_{n–1}\right)+(p_n–{\mu }_n{)}^2\\ & =& \left(p_n–{\mu }_n\right)\left(p_n–{\mu }_{n–1}\right)\end{array}$$and so

$$\begin{array}{ccc}{\mu }_n& =& {\mu }_{n–1}+\frac{1}{n}\left(p_n–{\mu }_{n–1}\right)\\ M_n& =& M_{n–1}+\left(p_n–{\mu }_n\right)\left(p_n–{\mu }_{n–1}\right)\\ {\sigma }_n^2& =& \frac{1}{n}{M}_n\end{array}$$This approach is known as Welford’s online algorithm and is numerically stable, O(1), and compuatble in one pass.

Biasedness of Estimation.

Under weak statistical assumption that the signal is a sequence of random variables with a stationary first moment $\mu$, then $E{\mu }_n=\mu$, that is the estimator unbiased. If the random variables are uncorrelated and have a stationary finite variance, then $E{\sigma }_n^2=\frac{n–1}{n}{\sigma }^2$, which is biased, but an easy fix can be made to obtain an unbiased estimator: $s_n^2=\frac{n}{n–1}{\sigma }_n^2$. This is typically not a serious issue since the estimators are all asymptotically unbiased. To see this, note that $$\begin{array}{ccc}E{M}_n& =& E\sum _{i=\mathrm{1..}n}{\left(p_i–{\mu }_n\right)}^2\\ & =& E\sum _{i=\mathrm{1..}n}{\left(p_i–\mu +\mu –{\mu }_n\right)}^2\\ & =& E\sum _{i=\mathrm{1..}n}\left[{\left(p_i–\mu \right)}^2+2\left(p_i–\mu \right)\left(\mu –{\mu }_n\right)+(\mu –{\mu }_n{)}^2\right]\\ & =& n{\sigma }^2–\frac{2}{n}\sum _{i=\mathrm{1..}n}\sum _{j–\mathrm{1..}n}E\left(p_i–\mu \right)\left(p_j–\mu \right)+{\sigma }^2\\ & =& n{\sigma }^2–\frac{2}{n}n{\sigma }^2+{\sigma }^2\\ & =& \left(n–1\right){\sigma }^2\end{array}$$