Online fitting for normal distributions

Question

I was wondering if there exist efficient online or dynamic algorithm for fitting a normal distribution to data as it comes in. I am interested in two variants:

1. The algorithm is fed data points one at a time, and has to update its previous best-fit to account for the new point at each step.

2. The algorithm is fed $n$ data points initially. At each time step one data point is removed (the oldest one of the $n$) and one data point is added. The algorithm has to efficiently update the best-fit.

Bonus points if you know an implementation of this algorithm in Matlab. Note, that the default Matlab normfit is not online, and it would be too computationally intensive to have to refit the data with it at every timestep.

Answer 1

For both variants you need to compute the variance and mean of your data to estimate the parameters of your normal distribution.

For number 1 Wikipedia cites a stable algorithm (and included pseudocode!) from The Art of Computer Programming, volume 2: Seminumerical Algorithms, 3rd edition that returns both here.

Number 2 is slightly trickier to program, but still straightforward. You will need a queue data structure to hold the incoming values. As each value comes in you dequeue the nth oldest value and queue the newest one. The mean is updated by subtracting the dequeued value from the sum and adding the new one. The variance is similarly calculated, but you need to do this with the sum of squares too. Just make sure the algorithm is stable before you have n numbers, and keep in mind that large n values may cause overflow issues when computing the sum of squares.

Answer 2

Here is the quick and dirty version, based on the sum of squares method, which is numerically inferior to Welford's method, but a lot faster to implement. I make use of an (anonymous) online sum function which takes vector $x$ and computes a vector $s$ such that $s_i = \sum_{\max(1,i-b+1) \le j \le i} x_j$.

A 'real' implementation of this function would be properly vectorized (i.e. would accept multidimensional array input and optional dimension argument, in analogy to Matlab's sum or cumsum functions), would use Welford's method, would do the right thing for nans, etc.

function [mu,sg] = online_normfit(x,boxwin);
%% online_normfit: fit mean and stdev of X 'online'
% 
% [mu,sg] = online_normfit(x,boxwin);
%
% takes an n vector x, and returns n vectors mu and sg
% such that
% mu(i) = mean(x(max(1,i-boxwin+1):i))
% sg(i) = std(x(max(1,i-boxwin+1):i))
%
% boxwin defaults to the length of x;
%
% nb: assumes there are no nans in x! 
% nb: this code can uderflow on variance computation because it uses
% the difference of squares method, not Welford's method. c.f.
% http://www.johndcook.com/blog/2008/09/26/comparing-three-methods-of-computing-standard-deviation/

if (~exist('boxwin','var')) || isempty(boxwin)
  boxwin = numel(x);
end

f_boxdif = @(cz,w)([cz(1:w);cz((w+1):end) - cz(1:(end-w))]);
f_boxsum = @(z,w)(f_boxdif(cumsum(z),w));

x1 = f_boxsum(x,boxwin);   % sums of x
x2 = f_boxsum(x.^2,boxwin);  % sums of x^2
nn = min((1:numel(x1))',boxwin);  % number of observations in each sum

mu = x1 ./ nn;
% this can underflow! also the first value is inf!
vr = (x2 - nn .* (mu.^2)) ./ (nn - 1);
sg = sqrt(vr);

end %function