Numerically stable "cumulative variance" function in O(n)

[stevengj]

@maxkapur

But (if I am understanding this code correctly) I think that this code ultimately relies on the computation mean(X^2) - mean(X)^2, and that computation itself is numerically unstable even if the two terms are calculated to a high degree of accuracy.

Why not use something like Welford’s algorithm?

function cumvar(x)
    sumsqr = (prevmean = zero(eltype(x)) / 1)^2
    v = Vector{typeof(sumsqr)}(undef, length(x))
    for (n, x) in enumerate(x)
        mean = (prevmean*(n-1) + x) / n
        sumsqr += (x - prevmean) * (x - mean)
        v[n] = sumsqr / n
        prevmean = mean
    end
    return v
end