using Printf

using Distributions
using StatsPlots
default(fmt=:png, size=(500, 300),
    titlefontsize=10, tickfontsize=6, guidefontsize=9,
    plot_titlefontsize=10)

# Sterneの信頼区間の定義で使用
using Roots
using StatsFuns: logit

# highest density interval (HDI) の構成で使用
using Optim

plot(x -> cdf(Binomial(10, 0.3), x), -1.5, 11.5; label="")
plot!(xtick=-1:11, ytick=0:0.05:1)
plot!(xguide="x", yguide="p")
title!("p = cdf(Binomial(10, 0.3), x)")

P1 = plot(z -> cdf(Normal(0,1), z), -4.5, 4.5; label="")
plot!(xtick=-10:10, ytick=0:0.1:1)
plot!(xguide="z", yguide="p")
title!("p = cdf(Normal(0,1), z)")

P2 = plot(p -> quantile(Normal(0,1), p), 0, 1; label="")
plot!(ylim=(-4.5, 4.5))
plot!(ytick=-10:10, xtick=0:0.1:1)
plot!(yguide="z", xguide="p")
title!("z = quantile(Normal(0,1), p)")

plot(P1, P2; size=(640, 320))

x = range(-4.5, 4.5, 1000)
plot(p -> quantile(Normal(0,1), p), 0, 1; label="quantile")
plot!(cdf.(Normal(0,1), x), x; label="inverse of cdf", ls=:dash)
plot!(ylim=(-4.5, 4.5))
plot!(ytick=-10:10, xtick=0:0.1:1)
plot!(yguide="z", xguide="p")
title!("Normal(0,1)")
plot!(size=(320, 320), legend=:topleft)

P1 = plot(z -> cdf(Binomial(10, 0.3), z), -0.5, 10.5; label="")
plot!(xtick=-1:100, ytick=0:0.1:1)
plot!(xguide="z", yguide="p")
title!("p = cdf(Binomial(10, 0.3), z)")

P2 = plot(p -> quantile(Binomial(10, 0.3), p), 0, 1; label="")
plot!(ylim=(-0.5, 10.5))
plot!(ytick=-1:100, xtick=0:0.1:1)
plot!(yguide="z", xguide="p")
title!("z = quantile(Binomial(10, 0.3), p)")

plot(P1, P2; size=(640, 320))

x = range(-0.5, 10.5, 1000)
plot(p -> quantile(Binomial(10, 0.3), p), 0, 1; label="quantile")
plot!(cdf.(Binomial(10, 0.3), x), x; label="inverse of cdf", ls=:dash)
plot!(ylim=(-0.5, 10.5))
plot!(ytick=-1:100, xtick=0:0.1:1)
plot!(yguide="z", xguide="p")
title!("Binomial(10, 0.3)")
plot!(size=(320, 320), legend=:topleft)

@show c = cquantile(Normal(0,1), 0.05/2)
plot(Normal(0,1), -4.5, 4.5; label="", xtick=-5:5, xguide="x")
plot!(Normal(0,1), c, 4.5;
    fillrange=0, c=1, fc=:pink, label="0.05/2 = 2.5%")
vline!([c]; label="x = $(round(c; digits=2))", c=:red, ls=:dash)
title!("Normal(0,1)"; legend=:topleft)

@show c = cquantile(Normal(0,1), 0.05/2)
plot(Normal(0,1), -4.5, 4.5; label="", xtick=-5:5, xguide="x")
plot!(Normal(0,1), c, 4.5;
    fillrange=0, c=1, fc=:pink, label="5%")
plot!(Normal(0,1), -4.5, -c;
    fillrange=0, c=1, fc=:pink, label="")
vline!([c]; label="x = $(round(c; digits=2))", c=:red, ls=:dash)
title!("Normal(0,1)"; legend=:topleft)

safediv(x, y) = x == 0 ? x : y == Inf ? zero(y) : x/y

function pvalue_wilson(k, n, p)
    p̂ = k/n
    se = √(p*(1 - p)/n)
    2ccdf(Normal(0,1), safediv(abs(p̂ - p), se))
end

function confint_wilson(k, n, α)
    z = cquantile(Normal(0,1), α/2)
    p̂ = k/n
    a, b, c = 1 + z^2/n, p̂ + z^2/(2n), p̂^2
    sqrtD = √(b^2 - a*c)
    (b - sqrtD)/a, (b + sqrtD)/a
end

@show n, k, α = 10, 3, 0.05
@show ci = confint_wilson(k, n, α) |> collect

P_Wilson = plot(p -> pvalue_wilson(k, n, p), 0, 1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI")
title!("Wilson's P-value function for n=$n, k=$k")

n = 10
k = 0:n
p = range(0, 1, 200)
heatmap(k, p, (k,p)->pvalue_wilson(k, n, p); clim=(0,1))
plot!(xtick=0:n, ytick=0:0.1:1)
plot!(xguide="k", yguide="p")
plot!(size=(350, 320))
title!("P-value (n=$n)")

P_Wilson

function pvalue_wald(k, n, p)
    p̂ = k/n
    sehat = √(p̂*(1 - p̂)/n)
    2ccdf(Normal(0,1), abs(p̂ - p)/sehat)
end

function confint_wald(k, n, α)
    z = cquantile(Normal(0,1), α/2)
    p̂ = k/n
    sehat = √(p̂*(1 - p̂)/n)
    p̂ - z*sehat, p̂ + z*sehat
end

@show n, k, α = 10, 3, 0.05
@show ci = confint_wald(k, n, α) |> collect

P_Wald = plot(p -> pvalue_wald(k, n, p), -0.1, 1.1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI")
title!("Wald's P-value function for n=$n, k=$k")

mypdf(dist::DiscreteUnivariateDistribution, x) = pdf(dist, round(Int, x))

function plot_bin_normal(n, p, p̂;
        μ = n*p, σ = √(n*p*(1-p)), σ̂ = √(n*p̂*(1-p̂)),
        xlim = (max(-1, μ - 4.5σ), min(n+1, μ + 4.5σ)),
        titlehead = "", kwargs...)
    bin = Binomial(n, p)
    normal = Normal(μ, σ̂)
    title = titlehead *
        "Bin(n=$n, p=$(@sprintf("%.2f", p))), p̂=$(@sprintf("%.2f", p̂))"
    plot(; legend=:outertop, legendfontsize=10, title)
    plot!(x -> mypdf(bin, x), xlim...; label="")
    plot!(normal, xlim...; label="", ls=:dash)
    plot!(; kwargs...)
end

n, p̂ = 30, 0.3
PP = [plot_bin_normal(n, p, p̂; xlim=(-1, n+1), titlefontsize=8)
    for p in 0.1:0.1:0.9]
plot(PP...; layout=(3, 3), size=(700, 500))

n, p̂ = 30, 0.3
ps = 0.05:0.01:0.95
ps = [fill(ps[begin], 10); ps; fill(ps[end], 10); reverse(ps)]
anim = @animate for p in ps
    plot_bin_normal(n, p, p̂;
        xlim=(-1, n+1), ylim=(-0.01, 0.35),
        xtick=0:n, ytick=0:0.05:1,
        titlehead="Wald: ")
end
gif(anim, "Wald_normal_approx.gif")

n = 30
ps = 0.05:0.01:0.95
ps = [fill(ps[begin], 10); ps; fill(ps[end], 10); reverse(ps)]
anim = @animate for p in ps
    plot_bin_normal(n, p, p;
        xlim=(-1, n+1), ylim=(-0.01, 0.35),
        xtick=0:n, ytick=0:0.05:1,
        titlehead="Wilson: ")
end
gif(anim, "Wilson_normal_approx.gif")

n, k = 10, 3
plot(p -> ccdf(Binomial(n,p), k-1);
    label="ccdf(Binomial(n,p), k-1)")
plot!(p -> cdf(Beta(k, n-k+1), p), 0, 1;
    label="cdf(Beta(k, n-k+1), p)", ls=:dash)
plot!(legend=:bottomright)
plot!(xguide="p")
plot!(xtick=0:0.1:1, ytick=0:0.1:1)
title!("n=$n, k=$k")

n, k = 10, 3
plot(p -> cdf(Binomial(n,p), k);
    label="cdf(Binomial(n,p), k)")
plot!(p -> ccdf(Beta(k+1, n-k), p), 0, 1;
    label="ccdf(Beta(k+1, n-k), p)", ls=:dash)
plot!(legend=:topright)
plot!(xguide="p")
plot!(xtick=0:0.1:1, ytick=0:0.1:1)
title!("n=$n, k=$k")

function pvalue_cp(k, n, p)
    bin = Binomial(n, p)
    min(1, 2cdf(bin, k), 2ccdf(bin, k-1))
end

function confint_cp(k, n, α)
    quantile(Beta(k, n-k+1), α/2), cquantile(Beta(k+1, n-k), α/2)
end

@show n, k, α = 10, 3, 0.05
@show ci = confint_cp(k, n, α) |> collect

P_CP = plot(p -> pvalue_cp(k, n, p), 0, 1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI")
title!("Clopper-Pearson P-value function for n=$n, k=$k")

x ⪅ y = x < y || x ≈ y

function pvalue_sterne_naive(k, n, p)
    bin = Binomial(n, p)
    sum(pdf(bin, i) for i in support(bin) if pdf(bin, i) ⪅ pdf(bin, k))
end

@time pvalue_sterne_naive(3*10^6, 10*10^6, 0.301)
@time pvalue_sterne_naive(3*10^6, 10*10^6, 0.301)
@time pvalue_sterne_naive(3*10^6, 10*10^6, 0.301)

x ⪅ y = x < y || x ≈ y

_pdf_le(x, (dist, y)) =  pdf(dist, x) ⪅ y

function _search_boundary(f, x0, Δx, param)
    x = x0
    if f(x, param)
        while f(x - Δx, param) x -= Δx end
    else
        x += Δx
        while !f(x, param) x += Δx end
    end
    x
end

function pvalue_sterne(dist::DiscreteUnivariateDistribution, x)
    Px = pdf(dist, x)
    Px == 0 && return Px
    Px == 1 && return Px
    m = mode(dist)
    Px ≈ pdf(dist, m) && return one(Px)
    if x < m
        y = _search_boundary(_pdf_le, 2m - x, 1, (dist, Px))
        cdf(dist, x) + ccdf(dist, y-1)
    else # x > m
        y = _search_boundary(_pdf_le, 2m - x, -1, (dist, Px))
        cdf(dist, y) + ccdf(dist, x-1)
    end
end

function pvalue_sterne(k, n, p)
    pvalue_sterne(Binomial(n, p), k)
end

function confint_sterne(k, n, α)
    a, b = confint_cp(k, n, α/10)
    ps = find_zeros(a-√eps(), b+√eps()) do p
        logit(0 < p ≤ 1 ? pvalue_sterne(k, n, p) : zero(p)) - logit(α)
    end
    first(ps), last(ps)
end

@time pvalue_sterne(3*10^6, 10*10^6, 0.301)
@time pvalue_sterne(3*10^6, 10*10^6, 0.301)
@time pvalue_sterne(3*10^6, 10*10^6, 0.301)

@time confint_sterne(3*10^6, 10*10^6, 0.001)
@time confint_sterne(3*10^6, 10*10^6, 0.001)
@time confint_sterne(3*10^6, 10*10^6, 0.001)

@show n, k, α = 10, 3, 0.05
@show ci = confint_sterne(k, n, α) |> collect

P_Sterne = plot(p -> pvalue_sterne(k, n, p), 0, 1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI")
title!("Sterne's P-value function for n=$n, k=$k")

plot(P_Wilson, P_Wald, P_CP, P_Sterne; layout=(2, 2), size=(800, 500))

@show n, k, α = 10, 3, 0.05
@show confint_wilson(k, n, α)
@show confint_wald(k, n, α)
@show confint_cp(k, n, α)
@show confint_sterne(k, n, α)

plot()
plot!(p -> pvalue_wilson(k, n, p), -0.1, 0.8; label="Wilson")
plot!(p -> pvalue_wald(k, n, p); label="Wald", ls=:dash)
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dashdot)
plot!(p -> pvalue_sterne(k, n, p); label="Sterne", ls=:dashdotdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 30, 9, 0.05
@show confint_wilson(k, n, α)
@show confint_wald(k, n, α)
@show confint_cp(k, n, α)
@show confint_sterne(k, n, α)

plot()
plot!(p -> pvalue_wilson(k, n, p), 0.05, 0.6; label="Wilson")
plot!(p -> pvalue_wald(k, n, p); label="Wald", ls=:dash)
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dashdot)
plot!(p -> pvalue_sterne(k, n, p); label="Sterne", ls=:dashdotdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 100, 30, 0.05
@show confint_wilson(k, n, α)
@show confint_wald(k, n, α)
@show confint_cp(k, n, α)
@show confint_sterne(k, n, α)

plot()
plot!(p -> pvalue_wilson(k, n, p), 0.15, 0.5; label="Wilson")
plot!(p -> pvalue_wald(k, n, p); label="Wald", ls=:dash)
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dashdot)
plot!(p -> pvalue_sterne(k, n, p); label="Sterne", ls=:dashdotdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 300, 90, 0.05
@show confint_wilson(k, n, α)
@show confint_wald(k, n, α)
@show confint_cp(k, n, α)
@show confint_sterne(k, n, α)

plot()
plot!(p -> pvalue_wilson(k, n, p), 0.2, 0.45; label="Wilson")
plot!(p -> pvalue_wald(k, n, p); label="Wald", ls=:dash)
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dashdot)
plot!(p -> pvalue_sterne(k, n, p); label="Sterne", ls=:dashdotdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

function hdi(dist::ContinuousUnivariateDistribution, α = 0.05; alg = Brent())
    f(t) = quantile(dist, t + (1 - α)) - quantile(dist, t)
    o = optimize(f, 0, α, alg)
    p = o.minimizer
    quantile.(dist, (p, p + (1 - α)))
end

function pvalue_hdi(dist::ContinuousUnivariateDistribution, x₀; xlim = extrema(dist))
    p₀ = pdf(dist, x₀)
    m = mode(dist)
    f(x) = pdf(dist, x) - p₀
    if x₀ == m
        1.0
    elseif x₀ > m
        x₁ = find_zero(f, (xlim[begin], m))
        cdf(dist, x₁) + ccdf(dist, x₀)
    else
        x₁ = find_zero(f, (m, xlim[end]))
        cdf(dist, x₀) + ccdf(dist, x₁)
    end
end

function pvalue_hdi(k, n, p; a=1, b=1)
    posterior = Beta(k+a, n-k+b)
    if k+a ≤ 1
        return ccdf(posterior, p)
    elseif n-k+b ≤ 1
        return cdf(posterior, p)
    end
    pvalue_hdi(posterior, p)
end

function credint_hdi(k, n, α; a=1, b=1)
    posterior = Beta(k+a, n-k+b)
    hdi(posterior, α)
end

@show n, k, α = 10, 3, 0.05
@show a, b = 1, 1
@show ci = credint_hdi(k, n, α; a, b) |> collect

P_HDI = plot(p -> pvalue_hdi(k, n, p; a, b), 0, 1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI (HDI)")
title!("Bayesian P-value function for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 10, 3, 0.05
@show a, b = 1, 1
@show ci = credint_hdi(k, n, α; a, b) |> collect

f(p) = pdf(Beta(k+a, n-k+b), p)
Q_HDI = plot(f, 0, 1; label="")
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.1:1)
plot!(ci, fill(f(ci[end]), 2); label="$(100(1-α))% CI (HDI)")
title!("Posterior for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 10, 3, 0.05
@show a, b = 1, 1

plot(p -> pvalue_hdi(k, n, p; a, b), 0, 1; label="Bayesian (HDI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 30, 9, 0.05
@show a, b = 1, 1

plot(p -> pvalue_hdi(k, n, p; a, b), 0.05, 0.6; label="Bayesian (HDI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 100, 30, 0.05
@show a, b = 1, 1

plot(p -> pvalue_hdi(k, n, p; a, b), 0.15, 0.5; label="Bayesian (HDI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

function pvalue_eti(k, n, p; a=1//3, b=1//3)
    posterior = Beta(k+a, n-k+b)
    min(1, 2cdf(posterior, p), 2ccdf(posterior, p))
end

function credint_eti(k, n, α; a=1//3, b=1//3)
    posterior = Beta(k+a, n-k+b)
    quantile(posterior, α/2), cquantile(posterior, α/2)
end

@show n, k, α = 10, 3, 0.05
@show a, b = 1/2, 1/2
@show ci = credint_eti(k, n, α; a, b) |> collect

P_ETI = plot(p -> pvalue_eti(k, n, p; a, b), 0, 1; label="")
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
plot!(ci, fill(α,2); label="$(100(1-α))% CI (ETI)")
title!("Bayesian P-value function for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 10, 3, 0.05
@show a, b = 1/2, 1/2

plot(p -> pvalue_eti(k, n, p; a, b), 0, 1; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 30, 9, 0.05
@show a, b = 1/2, 1/2

plot(p -> pvalue_eti(k, n, p; a, b), 0.05, 0.6; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 100, 30, 0.05
@show a, b = 1/2, 1/2

plot(p -> pvalue_eti(k, n, p; a, b), 0.15, 0.5; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 10, 3, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0, 1; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 30, 9, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0.05, 0.6; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k, α = 100, 30, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0.15, 0.5; label="Bayesian (ETI)")
plot!(p -> pvalue_wilson(k, n, p); label="Wilson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k = 10, 3
@show a, b = 1//3, 1//3

P1 = plot(Beta(k+a, n-k+b), 0, 0.85; label="Beta(k+a, n-k+b)")
plot!(Beta(k, n-k+1), 0, 0.85; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.1:1)
title!("n=$n, k=$k, a=$a, b=$b")

P2 = plot(Beta(k+a, n-k+b), 0., 0.85; label="Beta(k+a, n-k+b)")
plot!(Beta(k+1, n-k), 0, 0.85; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.1:1)
title!("n=$n, k=$k, a=$a, b=$b")

plot(P1, P2; size=(800, 250), leftmargin=4Plots.mm, bottommargin=4Plots.mm)

@show n, k, α = 10, 3, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0, 0.85; label="Bayesian (ETI)")
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.1:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k = 30, 9
@show a, b = 1//3, 1//3

P1 = plot(Beta(k+a, n-k+b), 0.05, 0.65; label="Beta(k+a, n-k+b)")
plot!(Beta(k, n-k+1), 0.05, 0.65; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

P2 = plot(Beta(k+a, n-k+b), 0.05, 0.65; label="Beta(k+a, n-k+b)")
plot!(Beta(k+1, n-k), 0.05, 0.65; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

plot(P1, P2; size=(800, 250), leftmargin=4Plots.mm, bottommargin=4Plots.mm)

@show n, k, α = 30, 9, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0.05, 0.65; label="Bayesian (ETI)")
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k = 100, 30
@show a, b = 1//3, 1//3

P1 = plot(Beta(k+a, n-k+b), 0.15, 0.5; label="Beta(k+a, n-k+b)")
plot!(Beta(k, n-k+1), 0.15, 0.5; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

P2 = plot(Beta(k+a, n-k+b), 0.15, 0.5; label="Beta(k+a, n-k+b)")
plot!(Beta(k+1, n-k), 0.15, 0.5; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

plot(P1, P2; size=(800, 250), leftmargin=4Plots.mm, bottommargin=4Plots.mm)

@show n, k, α = 100, 30, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0.15, 0.5; label="Bayesian (ETI)")
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

@show n, k = 300, 90
@show a, b = 1//3, 1//3

P1 = plot(Beta(k+a, n-k+b), 0.2, 0.45; label="Beta(k+a, n-k+b)")
plot!(Beta(k, n-k+1), 0.2, 0.45; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

P2 = plot(Beta(k+a, n-k+b), 0.15, 0.45; label="Beta(k+a, n-k+b)")
plot!(Beta(k+1, n-k), 0.15, 0.5; label="Beta(k, n-k+1)", ls=:dash)
plot!(xguide="p", yguide="probability density")
plot!(xtick=0:0.05:1)
title!("n=$n, k=$k, a=$a, b=$b")

plot(P1, P2; size=(800, 250), leftmargin=4Plots.mm, bottommargin=4Plots.mm)

@show n, k, α = 300, 90, 0.05
@show a, b = 1//3, 1//3

plot(p -> pvalue_eti(k, n, p; a, b), 0.2, 0.45; label="Bayesian (ETI)")
plot!(p -> pvalue_cp(k, n, p); label="Clopper-Pearson", ls=:dash)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:0.05:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k, a=$a, b=$b")

x ⪆ y = x > y || x ≈ y

function coverprob(pvaluefunc, n, p, α)
    bin = Binomial(n, p)
    sum(support(bin)) do k
        (pvaluefunc(k, n, p) ⪆ α) * pdf(bin, k)
    end
end

function plot_coverprob(n, α;
        prior_hdi = (1, 1), prior_eti = (1//3, 1//3),
        ylim = (1-2α, 1), lw = 1, kwargs...)
    @show n, α
    p = range(0, 1, 2001)

    CP_Wilson = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob(pvalue_wilson, n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Wilson: n=$n, α=$α")

    CP_Wald = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob(pvalue_wald, n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Wald: n=$n, α=$α")

    CP_CP = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob(pvalue_cp, n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Clopper-Pearson: n=$n, α=$α")

    CP_Sterne = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob(pvalue_sterne, n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Sterne: n=$n, α=$α")

    a, b = prior_hdi
    CP_HDI = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob((k,n,p)->pvalue_hdi(k,n,p; a,b), n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Bayesian HDI: prior=($a, $b), n=$n, α=$α")

    a, b = prior_eti
    CP_ETI = plot(; ylim, xtick=0:0.1:1, ytick=0:0.01:1)
    plot!(p, p -> coverprob((k,n,p)->pvalue_eti(k,n,p; a,b), n, p, α); label="", lw)
    hline!([1-α]; label="", ls=:dot)
    plot!(xguide="p")
    title!("Bayesian ETI: prior=($a, $b), n=$n, α=$α")

    plot(CP_CP, CP_Sterne, CP_Wilson, CP_Wald, CP_HDI, CP_ETI;
        size=(800, 1000), layout=(3, 2))
    plot!(plot_title="coverage probabilities")
    plot!(; kwargs...)
end

plot_coverprob(10, 0.05; ylim=(0.88, 1))

plot_coverprob(30, 0.05)

plot_coverprob(100, 0.05; lw=0.7)

plot_coverprob(300, 0.05; lw=0.5)

x ⪅ y = x < y || x ≈ y

function cdfpval(pvaluefunc, n, p, α)
    bin = Binomial(n, p)
    sum(support(bin)) do k
        (pvaluefunc(k, n, p) ⪅ α) * pdf(bin, k)
    end
end

function plot_cdfpval(n, p; f=Bool[1,1,1,1])
    α = range(0, 1, 1001)
    P1 = plot(legend=:topleft)
    f[1] && plot!(α, α->cdfpval(pvalue_wilson, n, p, α); label="Wilson", c=1)
    f[2] && plot!(α, α->cdfpval(pvalue_wald, n, p, α); label="Wald", ls=:dash, c=2)
    f[3] && plot!(α, α->cdfpval(pvalue_cp, n, p, α); label="Clopper-Pearson", ls=:dot, lw=1.3, c=3)
    f[4] && plot!(α, α->cdfpval(pvalue_sterne, n, p, α); label="Sterne", ls=:dashdot, lw=1.3, c=4)
    plot!([0,1], [0,1]; label="", c=:black, ls=:dot, lw=0.5)
    plot!(xtick=0:0.1:1, ytick=0:0.1:1)
    plot!(xguide="α", yguide="probabilty of P-value ≤ α")

    α = range(0, 0.1, 1001)
    P2 = plot(legend=:topleft)
    f[1] && plot!(α, α->cdfpval(pvalue_wilson, n, p, α); label="Wilson", c=1)
    f[2] && plot!(α, α->cdfpval(pvalue_wald, n, p, α); label="Wald", ls=:dash, c=2)
    f[3] && plot!(α, α->cdfpval(pvalue_cp, n, p, α); label="Clopper-Pearson", ls=:dot, lw=1.3, c=3)
    f[4] && plot!(α, α->cdfpval(pvalue_sterne, n, p, α); label="Sterne", ls=:dashdot, lw=1.3, c=4)
    plot!([0,0.1], [0,0.1]; label="", c=:black, ls=:dot, lw=0.5)
    plot!(xtick=0:0.01:1, ytick=0:0.01:1)
    plot!(xguide="α", yguide="probabilty of P-value ≤ α")

    plot!(P1, P2; size=(640, 320))
    plot!(plot_title="Assumption: data is generated by Binomial($n, $p).")
end

n, p = 10, 0.3
plot_cdfpval(n, p; f=Bool[1,1,0,0])

plot_cdfpval(n, p; f=Bool[0,0,1,1])

n, p = 30, 0.3
plot_cdfpval(n, p; f=Bool[1,1,0,0])

plot_cdfpval(n, p; f=Bool[0,0,1,1])

n, p = 100, 0.3
plot_cdfpval(n, p; f=Bool[1,1,0,0])

plot_cdfpval(n, p; f=Bool[0,0,1,1])

n, p = 300, 0.3
plot_cdfpval(n, p; f=Bool[1,1,0,0])

plot_cdfpval(n, p; f=Bool[0,0,1,1])

@show n, k, α = 10, 0, 0.05
@show a_eti, b_eti = 1//3, 1//3
@show a_hdi, b_hdi = 1, 1
pmin, pmax = 0, 0.56
pstep = 0.05

plot(p -> pvalue_wilson(k, n, p), pmin, pmax;
    label="Wilson score", ls=:solid)
plot!(p -> pvalue_cp(k, n, p), pmin, pmax;
    label="Clopper-Pearson", ls=:dot)
plot!(p -> pvalue_sterne(k, n, p), pmin, pmax;
    label="Sterne", ls=:dashdotdot)
plot!(p -> pvalue_hdi(k, n, p; a=a_hdi, b=b_hdi), pmin, pmax;
    label="Bayesian (HDI), prior=Beta($a_hdi, $b_hdi)", ls=:dash)
plot!(p -> pvalue_eti(k, n, p; a=a_eti, b=b_eti), pmin, pmax;
    label="Bayesian (ETI), prior=Beta($a_eti, $b_eti)", ls=:dashdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:pstep:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 40, 0, 0.05
@show a_eti, b_eti = 1//3, 1//3
@show a_hdi, b_hdi = 1, 1
pmin, pmax = 0, 0.16
pstep = 0.01

plot(p -> pvalue_wilson(k, n, p), pmin, pmax;
    label="Wilson score", ls=:solid)
plot!(p -> pvalue_cp(k, n, p), pmin, pmax;
    label="Clopper-Pearson", ls=:dot)
plot!(p -> pvalue_sterne(k, n, p), pmin, pmax;
    label="Sterne", ls=:dashdotdot)
plot!(p -> pvalue_hdi(k, n, p; a=a_hdi, b=b_hdi), pmin, pmax;
    label="Bayesian (HDI), prior=Beta($a_hdi, $b_hdi)", ls=:dash)
plot!(p -> pvalue_eti(k, n, p; a=a_eti, b=b_eti), pmin, pmax;
    label="Bayesian (ETI), prior=Beta($a_eti, $b_eti)", ls=:dashdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:pstep:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 160, 0, 0.05
@show a_eti, b_eti = 1//3, 1//3
@show a_hdi, b_hdi = 1, 1
pmin, pmax = 0, 0.04
pstep = 0.005

plot(p -> pvalue_wilson(k, n, p), pmin, pmax;
    label="Wilson score", ls=:solid)
plot!(p -> pvalue_cp(k, n, p), pmin, pmax;
    label="Clopper-Pearson", ls=:dot)
plot!(p -> pvalue_sterne(k, n, p), pmin, pmax;
    label="Sterne", ls=:dashdotdot)
plot!(p -> pvalue_hdi(k, n, p; a=a_hdi, b=b_hdi), pmin, pmax;
    label="Bayesian (HDI), prior=Beta($a_hdi, $b_hdi)", ls=:dash)
plot!(p -> pvalue_eti(k, n, p; a=a_eti, b=b_eti), pmin, pmax;
    label="Bayesian (ETI), prior=Beta($a_eti, $b_eti)", ls=:dashdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:pstep:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")

@show n, k, α = 640, 0, 0.05
@show a_eti, b_eti = 1//3, 1//3
@show a_hdi, b_hdi = 1, 1
pmin, pmax = 0, 0.01
pstep = 0.001

plot(p -> pvalue_wilson(k, n, p), pmin, pmax;
    label="Wilson score", ls=:solid)
plot!(p -> pvalue_cp(k, n, p), pmin, pmax;
    label="Clopper-Pearson", ls=:dot)
plot!(p -> pvalue_sterne(k, n, p), pmin, pmax;
    label="Sterne", ls=:dashdotdot)
plot!(p -> pvalue_hdi(k, n, p; a=a_hdi, b=b_hdi), pmin, pmax;
    label="Bayesian (HDI), prior=Beta($a_hdi, $b_hdi)", ls=:dash)
plot!(p -> pvalue_eti(k, n, p; a=a_eti, b=b_eti), pmin, pmax;
    label="Bayesian (ETI), prior=Beta($a_eti, $b_eti)", ls=:dashdot)
plot!(xguide="p", yguide="P-value")
plot!(xtick=0:pstep:1, ytick=0:0.05:1)
title!("P-value functions for n=$n, k=$k")