Monte Carlo Pi Error

pi_error.jl

function run_experiment(N; saverate=1)
    hits = 0
    dist = Uniform(-1,1)
    is = Int64[]
    pis = Float64[]
    
    for i in 1:N
        x, y = rand(dist), rand(dist)
        # speed-up computation by omiting sqrt
        if (x^2 + y^2) <= 1
            hits += 1
        end
        if (i-1)%saverate == 0
            push!(is, i)
            push!(pis, 4.0*hits/i)
        end
    end
    return is,pis
end

function plot_ensemble( ; Nrep=10, N=5000, saverate=1)
    f = Figure(resolution = (600,600), transparency=true)
    ax = Axis(f[1, 1], title="Quality of Approximation for $(Nrep) repetitions with $(Nrep) iterations")
    ylims!(ax, π-0.1, π+0.1)
    hlines!([π], color=:red, linewidth=1)
    
    for i=1:Nrep
        is, pis = run_experiment(N; saverate)
        lines!(is, pis, linewidth=1)
    end
    
    # statistical error analysis:
    # probability of hitting the circle is π/4
    # binomial distribution has variance: var(k) = n*p(1-p), where k = hitting the circle
    # hence, σ = sqrt(var(4*k/N)) = sqrt(4^2*var(k)/N^2) = 4*sqrt(N*p(1-p)/N^2) = 4*sqrt(p*(1-p)/N)
    p = pi/4
    σ(x) = 4*sqrt(p*(1-p)/x)
    
    # plot statistical error "envelopes"
    xs = 1:N
    ys = σ.(xs)
    # 95%
    lines!(xs, π .+ (2 .* ys), color=:red, linewidth=2)
    lines!(xs, π .- (2 .* ys), color=:red, linewidth=2)
    # 99%
    lines!(xs, π .+ (3 .* ys), color=:blue, linewidth=2)
    lines!(xs, π .- (3 .* ys), color=:blue, linewidth=2)

    save("pi-ensemble-plot.pdf",f)
    return f
end

Random.seed!(1234)
plot_ensemble(; Nrep=500, N=1_000_000, saverate=10_000)