Cov Float32 precision issues

[milankl]

Julia v1.6 on both mac and linux

julia> using Statistics
julia> A = 20*randn(Float32,10_000_000) .+ 100;
julia> cov(A,A)
293.97824f0

julia> var(A)
399.9234f0

While the true variance should be 400 for A the estimate from cov is way off. In contrast, var seems to be written in a way that the rounding error does not scale with the number of elements in A, as it should be.

[nalimilan]

Good catch. See https://github.com/JuliaLang/Statistics.jl/pull/85.

BTW, note that only the var method for arrays is accurate. If we use the method defined for general iterators, which does a single pass over the data (using Welford's algorithm), we get the same result as cov. Not sure something can be done about it.

julia> using Statistics

julia> A = 20*randn(Float32,10_000_000) .+ 100;

julia> cov(A,A)
394.65015f0

julia> var(A)
400.24118f0

julia> var((x for x in A))
394.6492f0

julia> var((x for x in A), mean=mean(A))
394.65015f0

[milankl]

Amazing work, thanks @nalimilan!