using Distributions, HypothesisTests, StatsBase, StatsPlots
theme(:solarized)


function calc_coverage(n, p, nrep; method = :wald)
    x   = rand(Binomial(n, p), nrep);
    bt  = BinomialTest.(x, n);
    CIs = confint.(bt; level = 0.95, tail = :both, method = method);

    res = Bool[]
    for i in eachindex(CIs)
        res = push!(res, CIs[i][1] .<= p && CIs[i][2] .>= p)
    end

    return(mean(res))
end


p = 0.2;
n = Int(1e4) .+ collect(0:100);
nrep = Int(1e7);

cov = zeros(length(n));
@time Threads.@threads for i in eachindex(n)
    cov[i] = calc_coverage(n[i], p, nrep; method = :clopper_pearson)
    if n[i] % 10 == 0
        println([n[i], cov[i]])
    end
end

plot(n, cov, label = "p = 0.2", xlabel = "n", ylabel = " Coverage")