Deleting Values in Welford’s Algorithm for Online Mean and Variance

To take a Gibbs sample, the previous sampled value for a variable is deleted from the sufficient statistics of the estimators that depend on it, the variable is resampled, and the the variable is added back into the estimators. What if our variable’s normal? Well, it turns out we can generalize Welford’s algorithm, which I describe in a comment to:

LingPipe Blog: Computing Sample Mean and Variance Online in One Pass

Here’s the basic algorithm, cribbed and simplified from John Cook’s site (see reference below), along with an extra method (unHandle(double)) added in for deletes:

private long mN = 0L;
private double mM = 0.0;
private double mS = 0.0;
public void handle(double x) {
     ++mN;
    double nextM = mM + (x – mM) / mN;
    mS += (x – mM) * (x – nextM);
    mM = nextM;
}
public void unHandle(double x) {
    if (mN == 0) {
        throw new IllegalStateException();
    } else if (mN == 1) { 
        mN = 0; mM = 0.0; mS = 0.0;
    } else {
        double mMOld = (mN * mM - x)/(mN-1);
        mS -= (x -  mM) * (x - mMOld);
        mM = mMOld;
        --mN;
    }
}
public double mean() {
    return mM;
}
public double varianceUnbiased() {
    return mN > 1 ? mS/(mN-1) : 0.0;
}

Works like a charm. I’m sure someone’s done this already, but I didn’t find any references in a quick check.

The same technique can obviously be applied to covariance and correlation calculations.

This’ll be in the next release of LingPipe.

References

Welford, B. P. 1962. Note on a method for calculating corrected sums of squares and products. Technometrics 4(3): 419-420.
Cook, John. Accurately computing running variance. Web page.

This entry was posted on July 7, 2009 at 11:32 am and is filed under Carp's Blog, LingPipe in Use, LingPipe News. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

11 Responses to “Deleting Values in Welford’s Algorithm for Online Mean and Variance”

Joseph Says:
August 9, 2010 at 5:30 am | Reply
I’m a student at University of Manitiba and had a question I thought you might be able to answer. I’m trying to identify and algorith that would be best suited to help me identify relationship between two arrays (each array generated froma separate archival story/article). Each array contains 1 column of keyword entries paired to a second column of frequencies that the keyword occurs within the original text.

The idea is that I should be able to compare these two arrays and identify if there are enough keywords in common to deem these two articles related.

Any suggestions on stats models or formulas to use? Point me in the right direction?

Thanks!
- lingpipe Says:
  August 9, 2010 at 12:32 pm | Reply
  This is not the place to ask a general question! We have a mailing list and e-mail (see the LingPipe home page).
  
  The usual thing to do is to treat the word to count maps as vectors with the words as the dimensions and then use standard vector cosine to compare them. This is all implemented in LingPipe, though has nothing to do with this post. Often, there’s a TF/IDF rescaling of the counts. Check out LingPipe’s class spell.TfIdfDistance for details and an implementation.
amocanu2013 Says:
November 18, 2014 at 4:59 pm | Reply
Maybe I’m missing sth, but according to http://www.johndcook.com/blog/standard_deviation/ you should have a case then mean == 1 in your handle() method. The part of the code I’m referring to is marked with this comment at that link // See Knuth TAOCP vol 2, 3rd edition, page 232
- Bob Carpenter Says:
  December 10, 2014 at 4:01 pm | Reply
  That’s just an efficiency thing, but I’m skeptical whether it would be faster or not. Branch mispredictions are costly, whereas not having the mN == 1 is a cost you only pay once per accumulation.
Steve Ash (@steve_ash) Says:
August 6, 2015 at 6:43 pm | Reply
I wish I had found this post before spending time coming up with my own implementation. I ran some simulations with mine and compared it to yours. The difference is negligible for anything that I care about, but if you’re curious, check out: http://blog.manycupsofcoffee.com/2015/08/streaming-standard-deviation-and.html
Matthew Fioravante Says:
August 20, 2015 at 11:48 am | Reply
You can easily extend this algorithm to compute correlation, covariance and all of the linear regression coefficients.

You just need to compute the cross product sum
add:
sxy += (x – mM) * (y – nextM);
remove:
sxy -= (x – mM) * (y – mMOld);

And then sxy / (n-1) is the covariance and so on…
- Matthew Fioravante Says:
  August 20, 2015 at 11:50 am | Reply
  Sorry types in my previous post:
  add:
  sxy += (x – mx) * (y – next_my);
  remove:
  sxy -= (x – mx) * (y – old_my);
  - Bob Carpenter Says:
    August 20, 2015 at 1:21 pm
    Although it’s C++, there’s a nice implementation of all of this in the Boost accumulators (though I don’t know if they give you the MLE divide-by-N estimate or the unbiased divide-by-(N-1) estimate. You can also accumulate Y^2 and use the definition of variance. It avoids all those potentially imprecise subtractions.
Marcus Nyberg Says:
April 27, 2017 at 3:48 am | Reply
Thanks! I tried doing this myself but ran into trouble. This works great!
biotinker Says:
March 24, 2020 at 1:16 pm | Reply
This code is incorrect and will not yield the mean of the input values.
Example in python:
>>> mN = 0
>>> mM = 0
>>> mS = 0
>>>
>>> def handle(x):
… global mN
… global mM
… global mS
… mN = mN + 1
… nextM = mM + (x – mM) / mN
… mS = mS + (x – mM) * (x – nextM)
… mM = nextM
… mN = mN + 1
… print(mM)
…
>>> handle(1)
1.0
>>> handle(1)
1.0
>>> handle(2)
1.2
>>> handle(3)
1.457142857142857
>>>
>>> realMean = (1 + 1 + 2 + 3)/4
>>> print(realMean)
1.75
- biotinker Says:
  March 24, 2020 at 1:29 pm | Reply
  Disregard the above, I have
  mN = mN + 1
  in there twice. Silly me