michael_frank_monte_carlo_study.jl

import Distributions: Binomial

function clt_ci(x)
    m = mean(x)
    s = std(x)
    n = length(x)
    se = s / sqrt(n)
    return max(0.0, m - 1.96 * se), min(1.0, m + 1.96 * se)
end

function bootstrap_ci(x, n_replicates = 1_000)
    means = Array(Float64, n_replicates)
    n = length(x)
    x_simulated = Array(Float64, n)
    for rep in 1:n_replicates
        for i in 1:n
            j = rand(1:n)
            x_simulated[i] = x[j]
        end
        means[rep] = mean(x_simulated)
    end
    lower = quantile(means, 0.025)
    upper = quantile(means, 0.975)
    return lower, upper
end

function simulate(n_sims, p)
    coverage_clt, coverage_bootstrap = 0, 0
    n_subs = 50
    n_obs = 10
    d = Array(Float64, n_subs)
    for sim in 1:n_sims
        for s in 1:n_subs
            d[s] = rand(Binomial(n_obs, p)) / n_obs
        end
        lower, upper = clt_ci(d)
        coverage_clt += Int(lower <= p <= upper)
        lower, upper = bootstrap_ci(d)
        coverage_bootstrap += Int(lower <= p <= upper)
    end
    return coverage_clt / n_sims, coverage_bootstrap / n_sims
end

julia> x, y = simulate(2_500, 0.90)
(0.9508,0.954)

julia> x, y = simulate(2_500, 0.99)
(0.864,0.9436)

julia> x, y = simulate(2_500, 0.999)
(0.3948,0.3936)