Testing summation error on distribution normalization

output.txt

size=1000
  Naive sum min: 0.9999999999999952 (error: 4.7739590058881731e-15)
  Naive sum max: 1.0000000000000056 (error: 5.5511151231257827e-15)
  KBN sum min  : 0.9999999999999999 (error: 1.1102230246251565e-16)
  KBN sum max  : 1.0000000000000000 (error: 0.0000000000000000e+00)
  Min p        : 0.0000000017882931 (       1.7882930846464359e-09)

size=10000
  Naive sum min: 0.9999999999999883 (error: 1.1657341758564144e-14)
  Naive sum max: 1.0000000000000142 (error: 1.4210854715202004e-14)
  KBN sum min  : 0.9999999999999999 (error: 1.1102230246251565e-16)
  KBN sum max  : 1.0000000000000000 (error: 0.0000000000000000e+00)
  Min p        : 0.0000000000376995 (       3.7699519733262038e-11)

size=100000
  Naive sum min: 0.9999999999999573 (error: 4.2743586448068527e-14)
  Naive sum max: 1.0000000000000273 (error: 2.7311486405778851e-14)
  KBN sum min  : 0.9999999999999999 (error: 1.1102230246251565e-16)
  KBN sum max  : 1.0000000000000000 (error: 0.0000000000000000e+00)
  Min p        : 0.0000000000002040 (       2.0404123948591026e-13)

size=1000000
  Naive sum min: 0.9999999999999227 (error: 7.7271522513910895e-14)
  Naive sum max: 1.0000000000000580 (error: 5.7953641885433171e-14)
  KBN sum min  : 0.9999999999999999 (error: 1.1102230246251565e-16)
  KBN sum max  : 1.0000000000000000 (error: 0.0000000000000000e+00)
  Min p        : 0.0000000000001159 (       1.1593101216101975e-13)

summation_error_test.jl

load("distributions.jl")
import Distributions

function naive_sum(a)
    s = 0
    for i = 1:length(a)
        s = s+a[i]
    end
    return s
end

dist = Exponential(1.5)
for sz in [1000 10000 100000 1000000]
    iters = min(1000, int(10000000/sz))
    a = zeros(iters);
    b = zeros(iters);
    min_sample = 1.0
    for i = 1:iters
        exp_samples = rand(dist, sz)
        sort!(>, exp_samples)
        exp_naive = exp_samples/naive_sum(exp_samples)
        exp_kbn = exp_samples/sum(exp_samples)
        if exp_naive[end] < min_sample
            min_sample = exp_naive[end]
        end
        if exp_kbn[end] < min_sample
            min_sample = exp_kbn[end]
        end
        a[i] = naive_sum(exp_naive)
        b[i] = sum(exp_kbn)
    end
    println()
    @printf("size=%d\n", sz)
    @printf("  Naive sum min: %1.16f (error: %1.16e)\n", min(a), 1.0-min(a))
    @printf("  Naive sum max: %1.16f (error: %1.16e)\n", max(a), max(a)-1.0)
    @printf("  KBN sum min  : %1.16f (error: %1.16e)\n", min(b), 1.0-min(b))
    @printf("  KBN sum max  : %1.16f (error: %1.16e)\n", max(b), max(b)-1.0)
    @printf("  Min p        : %1.16f (       %1.16e)\n", min_sample, min_sample)
end