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.

  • 3
    $\begingroup$ Online fitting of what? Just a marginal distribution? In that case, you would just use an online update of the mean and variance and these are trivial to obtain since $(\bar{X}_n, S^2_n)$ is Markov as a function of $n$. In fact, the main hurdle is to implement this in a way that avoids numerical cancellation. $\endgroup$ – cardinal Jan 31 '12 at 19:39
  • 2
    $\begingroup$ Normal distributions are usually fit by computing the mean and variance of the sample. Thus, the standard online methods to compute those statistics would work. $\endgroup$ – whuber Jan 31 '12 at 19:40
  • $\begingroup$ @whuber what should I do if I want to also compute the error on the mean, and error on the variance (as output by normfit) at some confidence level? $\endgroup$ – Artem Kaznatcheev Jan 31 '12 at 19:44
  • $\begingroup$ The standard error of the mean is computed from the variance and the SE of the variance is computed from the fourth central moment. The latter can be updated online in a manner similar to the other moments. In fact, once you can compute an online mean, you can use the same formula to compute an online mean square or mean fourth power, then combine those estimates in the usual way to obtain the central moments. $\endgroup$ – whuber Jan 31 '12 at 19:48
  • 2
    $\begingroup$ I guess @rm999 gave an answer already, so it would be inappropriate for me to expand the question. $\endgroup$ – Artem Kaznatcheev Jan 31 '12 at 19:55

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.

| cite | improve this answer | |

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);

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
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.