In [1]:
## Algorithms for the Variance

function badvar1(x::Vector{Float64})
    ex2 = 0.0
    ex = 0.0
    n = length(x)
    for i=1:n
      ex2 = ex2 + x[i]^2
      ex = ex + x[i]
    end
    @show ex2, ex^2
    return 1.0/(n-1)*(ex2 - (ex)^2/n)
end
Out[1]:
badvar1 (generic function with 1 method)
In [2]:
##
x = randn(100)
basevar = x -> (length(x)/(length(x)-1))*mean((x - mean(x)).^2)
@show badvar1(x)
@show basevar(x)
(ex2, ex ^ 2) = (114.37905866583934, 21.453908291731675)
badvar1(x) = 1.153176965484061
basevar(x) = 1.153176965484061
Out[2]:
1.153176965484061
In [3]:
##
x = randn(10000)+1e4
@show badvar1(x)
@show basevar(x)
(ex2, ex ^ 2) = (9.999998732727296e11, 9.999998632115106e15)
badvar1(x) = 1.006222521666229
basevar(x) = 1.0062227386805418
Out[3]:
1.0062227386805418
In [4]:
##
x = randn(10000)+1e6
@show badvar1(x)
@show basevar(x)
(ex2, ex ^ 2) = (9.99999987505691e15, 9.999999875046872e19)
badvar1(x) = 1.003900390039004
basevar(x) = 1.0072958627362085
Out[4]:
1.0072958627362085
In [5]:
##
x = randn(10000)+1e8
@show badvar1(x)
@show basevar(x)
(ex2, ex ^ 2) = (9.99999999842724e19, 9.999999998427243e23)
badvar1(x) = -3.2771277127712772
basevar(x) = 0.9839114279774787
Out[5]:
0.9839114279774787
In [6]:
##
function goodvar(x::Vector{Float64})
    n = length(x); mean = 0.0; m2 = 0.0; N = 0.0
    for i=1:n
        N = N + 1
        delta = x[i] - mean
        mean = mean + delta/N
        m2 = m2 + delta*(x[i]-mean)
    end
    return m2/(n-1)
end

x = randn(10000)+1e8
basevar = var
@show badvar1(x)
@show basevar(x)
@show goodvar(x)
(ex2, ex ^ 2) = (1.0000000002230891e20, 1.0000000002230964e24)
badvar1(x) = -72.0968096809681
basevar(x) = 0.9971540380862616
goodvar(x) = 0.9971540443441296