Implementing Welford's Algorithm (incremental variance calculation)

Hi. I'm looking to iteratively calculate variance since my home desktop doesn't have enough RAM. I've tried implementing the below algorithm (written in Python from Wikipedia) to generalize to n-dimension arrays (but I really only need n = 3), but I keep getting errors. Does anyone know of a Matlab implementation?

# for a new value newValue, compute the new count, new mean, the new M2.
# mean accumulates the mean of the entire dataset
# M2 aggregates the squared distance from the mean
# count aggregates the number of samples seen so far
def update(existingAggregate, newValue):
  (count, mean, M2) = existingAggregate
  count = count + 1 
  delta = newValue - mean
  mean = mean + delta / count
  delta2 = newValue - mean
  M2 = M2 + delta * delta2
return (count, mean, M2)

 # retrieve the mean, variance and sample variance from an aggregate
 def finalize(existingAggregate):
    (count, mean, M2) = existingAggregate
    (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1)) 
    if count < 2:
        return float('nan')
    else:
        return (mean, variance, sampleVariance)

Thanks!

[James Tursa]

Can you show us the MATLAB code you have so far?

Answer 1

[Jeff Miller]

RunStat on GitHub seems to have a MATLAB implementation (among others)

[Ersa U]

Awesome! The only issue is that it only accepts scalars while I need it to accept matricies.

Any idea how we could do that?

[Jeff Miller]

Not sure what you mean by "I need it to accept matricies".

If you get a whole batch of newValues at once (i.e., a matrix of them), then you can feed those into a RunStat accumulator one at a time using a for loop.

If you have many different kinds of newValues to be treated separately (i.e., each matrix position is a separate variable), then you can set up a separate accumulator for each matrix position and feed each accumulator its new value from each matrix (for loop again).

If you have many different kinds of newValues and want to accumulate some kind of variance/covariance matrix to reflect not only their individual variances but also their correlations, then the answer to your question is, "No, no idea, sorry." I don't know if there is a generalization of Welford's algorithm for accumulating covariances in a numerically stable way.