Element-wise variance of an iterator

Question

What's a numerically-stable way of taking the variance of an iterator elementwise? As an example, I would like to do something like

var((rand(4,2) for i in 1:10))

and get back a (4,2) matrix which is the variance in each coefficient. This throws an error using Julia's Base var. Is there a package that can handle this? Or an easy (and storage-efficient) way to do this using the Base Julia function? Or does one need to be developed on its own?

Answer 1

I went ahead and implemented a Welford algorithm to calculate this:

# Welford algorithm
# https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
function componentwise_meanvar(A;bessel=true)
  x0 = first(A)
  n = 0
  mean = zero(x0)
  M2 = zero(x0)
  delta = zero(x0)
  delta2 = zero(x0)
  for x in A
    n += 1
    delta .= x .- mean
    mean .+= delta./n
    delta2 .= x .- mean
    M2 .+= delta.*delta2
  end
  if n < 2
    return NaN
  else
    if bessel
      M2 .= M2 ./ (n .- 1)
    else
      M2 .= M2 ./ n
    end
    return mean,M2
  end
end

A few other algorithms are implemented in DiffEqMonteCarlo.jl as well. I'm surprised I couldn't find a library for this, but maybe will refactor this out someday.