Notebook

Bernoulli分布モデル¶

黒木玄

2020-02-11～2020-02-14

1 手計算のノート
2 共役事前分布
3 最尤法とベイズ法の予測分布
4 Kullback-Leibler情報量による予測誤差の定義
5 予測誤差の経験分布
6 平均予測誤差と平均二乗誤差
7 AICとWAICとLOOCV
8 BICとWBICと自由エネルギー
- 8.1 BICとWBICと自由エネルギーの関係
- 8.2 BICとWBICと自由エネルギーによって正しいモデル選択に失敗する確率
9 ベイズ統計と仮説検定や信頼区間との関係
- 9.1 P値函数およびその類似物の定義
- 9.2 P値函数のアニメーション

In [1]:

ENV["LINES"] = 200

using Distributions

using Plots
pyplot(fmt=:svg)
pyplotclf() = if backend() == Plots.PyPlotBackend(); PyPlot.clf(); end

using Base64
displayfile(mime, file; tag="img") = open(file) do f
    display("text/html", """<$tag src="data:$mime;base64,$(base64encode(f))"/>""")
end

using Printf
using Roots
using SpecialFunctions

rd(x, digits=3) = round(x; digits=digits)
safemul(x, y) = iszero(x) ? x : x*y
safexlogy(x, y) = iszero(x) ? x : x*log(y)

Out[1]:

safexlogy (generic function with 1 method)

手計算のノート¶

以下の手書きのノートで間違っていたとしても, コードの方では直っている可能性がある.

In [2]:

displayfile("image/png", "Bernoulli model.png")

共役事前分布¶

In [3]:

# mean(Beta(a, b)) → a/(a+b)
mean(Beta(3, 7))

Out[3]:

0.3

In [4]:

# lim_{b→0} mean(a, b) = a/a = 1
mean(Beta(1, 1e-8)), mean(Beta(2, 1e-8))

Out[4]:

(0.9999999900000002, 0.999999995)

In [5]:

# the improper "Beta(0, 0)" distribution
plot(x -> 1/(x*(1-x)), 0, 1; ylim=(-0.1, 100), label="1/(x(1-x))", size=(400, 250))

Out[5]:

In [6]:

# the improper "Beta(1, 0)" distribution
plot(x -> 1/(1-x), 0, 1; ylim=(-0.1, 100), label="1/(1-x)", size=(400, 250))

Out[6]:

In [7]:

# Jeffreys' prior
plot(x -> pdf(Beta(1/2, 1/2), x), 0, 1; ylim=(-0.1, 10), label="Beta(1/2, 1/2)", size=(400, 250))

Out[7]:

最尤法とベイズ法の予測分布¶

In [8]:

pred_MLE(n, k) = Bernoulli(k/n) # the predictive distribution of maximum likelihood estimation 
pred_Bayes(n, k; β=1.0, a=0.5, b=0.5) = Bernoulli((β*k+a)/(β*n+a+b)) # a=b=0.5 for the Jeffreys prior

n = 10
@show [pred_MLE(n, k).p   |> rd for k in 0:n]
@show [pred_Bayes(n, k).p |> rd for k in 0:n]
@show [pred_Bayes(n, k; a=1.0, b=1.0).p |> rd for k in 0:n]
@show [pred_Bayes(n, k; a=0.01, b=0.01).p |> rd for k in 0:n]
@show [pred_Bayes(n, k; β=0.0).p |> rd for k in 0:n]
@show [pred_Bayes(n, k; β=1e8).p |> rd for k in 0:n]
;

[(pred_MLE(n, k)).p |> rd for k = 0:n] = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
[(pred_Bayes(n, k)).p |> rd for k = 0:n] = [0.045, 0.136, 0.227, 0.318, 0.409, 0.5, 0.591, 0.682, 0.773, 0.864, 0.955]
[(pred_Bayes(n, k; a=1.0, b=1.0)).p |> rd for k = 0:n] = [0.083, 0.167, 0.25, 0.333, 0.417, 0.5, 0.583, 0.667, 0.75, 0.833, 0.917]
[(pred_Bayes(n, k; a=0.01, b=0.01)).p |> rd for k = 0:n] = [0.001, 0.101, 0.201, 0.3, 0.4, 0.5, 0.6, 0.7, 0.799, 0.899, 0.999]
[(pred_Bayes(n, k; β=0.0)).p |> rd for k = 0:n] = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
[(pred_Bayes(n, k; β=1.0e8)).p |> rd for k = 0:n] = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]

Kullback-Leibler情報量による予測誤差の定義¶

In [9]:

G(q, p) = -sum(safemul(pdf(q, x), logpdf(p, x)) for x in support(q)) # generalization error (or loss)
S(q) = G(q, q) # Shannon's information
KL(q, p) = G(q, p) - S(q) # Kullback-Leibler information, or prediction error of simulation of q by p

[(p, KL(Bernoulli(0.3), Bernoulli(p))) for p in 0:0.1:1]

Out[9]:

11-element Array{Tuple{Float64,Float64},1}:
 (0.0, Inf)                 
 (0.1, 0.15366358680379855) 
 (0.2, 0.028167557595283332)
 (0.3, 0.0)                 
 (0.4, 0.02160085414354651) 
 (0.5, 0.08228287850505178) 
 (0.6, 0.1837868973868122)  
 (0.7, 0.33891914415488134) 
 (0.8, 0.5826853020432395)  
 (0.9, 1.0325534177382862)  
 (1.0, Inf)

予測誤差の経験分布¶

In [10]:

make_samples(q::Bernoulli, n, L) = rand(Binomial(n, q.p), L)
ecdf(sample, x) = mean(sample .≤ x)

function plot_esitmation_errors(true_prob, n; L=10^6, bin=0:0.02:0.3, β=1.0, a=0.5, b=0.5)
    q = Bernoulli(true_prob)
    sample = make_samples(q, n, L)
    p_MLE   = pred_MLE.(n, sample)
    p_Bayes = pred_Bayes.(n, sample; β=β, a=a, b=b)
    KL_MLE   = KL.(q, p_MLE)
    KL_Bayes = KL.(q, p_Bayes)
    @show true_prob
    @show n
    
    P1 = plot(title="normalized histogram", titlefontsize=8)
    histogram!(KL_MLE;   bin=bin, label="KL_MLE",   alpha=0.5, norm=true)
    histogram!(KL_Bayes; bin=bin, label="KL_Bayes(β=$(β), a=$a, b=$b)", alpha=0.2, norm=true)
    x = range(0.01maximum(bin), maximum(bin), length=400)
    plot!(x, 2n*pdf.(Chisq(1), 2n*x);    label="2n pdf(2nx|χ²(1))", ls=:dot, color=:black)

    P2 = plot(title="empirical cumulative distribution function", titlefontsize=8)
    plot!(x, ecdf.(Ref(KL_MLE),   x); label="KL_MLE")
    plot!(x, ecdf.(Ref(KL_Bayes), x); label="KL_Bayes(β=$(β), a=$a, b=$b)", ls=:dash)
    x = range(extrema(bin)..., length=400)
    plot!(x, cdf.(Chisq(1), 2n*x);    label="cdf(2nx|χ²(1))", ls=:dot, color=:black)

    sleep(0.1)
    plot(P1, P2; size=(800, 250), layout=(1,2))
end

Out[10]:

plot_esitmation_errors (generic function with 1 method)

In [11]:

plot_esitmation_errors(0.3, 10)

true_prob = 0.3
n = 10

Out[11]:

In [12]:

plot_esitmation_errors(0.3, 100; bin=0:0.0025:0.05)

true_prob = 0.3
n = 100

Out[12]:

平均予測誤差と平均二乗誤差¶

In [13]:

MSE(q, p) = sum((p.p - q.p)^2*pdf(q, k) for k in support(q)) # mean squared error

# expextation value
E(dist, f) = sum(f(x)*pdf(dist, x) for x in support(dist))

# expextation value (extremum cases eliminated)
function EE(dist, f)
    s = support(dist)
    sum(f(x)*pdf(dist, x) for x in s[2:end-1])/(1 - pdf(dist, s[1]) - pdf(dist, s[end]))
end

function print_estimation_errors(true_prob, n; β=1.0, a=0.5, b=0.5)
    q = Bernoulli(true_prob)
    q_sample = Binomial(n, q.p)
    f_KL_Bayes(x) = KL(q, pred_Bayes(n, x; β=β, a=a, b=b))
    f_KL_MLE(x) = KL(q, pred_MLE(n, x))
    f_MSE_Bayes(x) = MSE(q, pred_Bayes(n, x; β=β, a=a, b=b))
    f_MSE_MLE(x) = MSE(q, pred_MLE(n, x))

    @show true_prob
    @show n
    @show β
    @show a
    @show b
    println()
    @show E(q_sample, f_KL_Bayes)
    @show E(q_sample, f_KL_MLE)
    @show EE(q_sample, f_KL_Bayes)
    @show EE(q_sample, f_KL_MLE)
    println()
    @show E(q_sample, f_MSE_Bayes)
    @show E(q_sample, f_MSE_MLE)
    @show EE(q_sample, f_MSE_Bayes)
    @show EE(q_sample, f_MSE_MLE)
    nothing
end

Out[13]:

print_estimation_errors (generic function with 1 method)

In [14]:

print_estimation_errors(0.3, 10)

true_prob = 0.3
n = 10
β = 1.0
a = 0.5
b = 0.5

E(q_sample, f_KL_Bayes) = 0.047718404765104
E(q_sample, f_KL_MLE) = Inf
EE(q_sample, f_KL_Bayes) = 0.03895091197067187
EE(q_sample, f_KL_MLE) = 0.05017959517455092

E(q_sample, f_MSE_Bayes) = 0.017685950413223135
E(q_sample, f_MSE_MLE) = 0.020999999999999994
EE(q_sample, f_MSE_Bayes) = 0.016314096917360407
EE(q_sample, f_MSE_MLE) = 0.018991401589615813

In [15]:

print_estimation_errors(0.3, 100)

true_prob = 0.3
n = 100
β = 1.0
a = 0.5
b = 0.5

E(q_sample, f_KL_Bayes) = 0.004979051813238021
E(q_sample, f_KL_MLE) = Inf
EE(q_sample, f_KL_Bayes) = 0.004979051813237705
EE(q_sample, f_KL_MLE) = 0.005112075811986031

E(q_sample, f_MSE_Bayes) = 0.0020625428879521616
E(q_sample, f_MSE_MLE) = 0.0020999999999999994
EE(q_sample, f_MSE_Bayes) = 0.0020625428879521347
EE(q_sample, f_MSE_MLE) = 0.0020999999999999712

In [16]:

function plot_mean_estimation_errors(n; 
        legend_KL=true, legend_MSE=true, 
        β₁=1.0, a₁=0.5, b₁=0.5,
        β₂=1.0, a₂=1.0, b₂=1.0
    )
    f_KL_MLE(true_prob, x)    = KL(Bernoulli(true_prob), pred_MLE(n, x))
    f_KL_BayesJ(true_prob, x) = KL(Bernoulli(true_prob), pred_Bayes(n, x; β=β₁, a=a₁, b=b₁))
    f_KL_Bayes(true_prob, x)  = KL(Bernoulli(true_prob), pred_Bayes(n, x; β=β₂, a=a₂, b=b₂))

    f_MSE_MLE(true_prob, x)    = MSE(Bernoulli(true_prob), pred_MLE(n, x))
    f_MSE_BayesJ(true_prob, x) = MSE(Bernoulli(true_prob), pred_Bayes(n, x; β=β₁, a=a₁, b=b₁))
    f_MSE_Bayes(true_prob, x)  = MSE(Bernoulli(true_prob), pred_Bayes(n, x; β=β₂, a=a₂, b=b₂))
    
    EE_KL_MLE(true_prob)    = EE(Binomial(n, true_prob), x->f_KL_MLE(true_prob, x))
    EE_KL_BayesJ(true_prob) = EE(Binomial(n, true_prob), x->f_KL_BayesJ(true_prob, x))
    EE_KL_Bayes(true_prob)  = EE(Binomial(n, true_prob), x->f_KL_Bayes(true_prob, x))

    E_MSE_MLE(true_prob)    = EE(Binomial(n, true_prob), x->f_MSE_MLE(true_prob, x))
    E_MSE_BayesJ(true_prob) = EE(Binomial(n, true_prob), x->f_MSE_BayesJ(true_prob, x))
    E_MSE_Bayes(true_prob)  = EE(Binomial(n, true_prob), x->f_MSE_Bayes(true_prob, x))

    true_prob = 0.01:0.01:0.99
    
    P1 = plot(title="mean prediction error (extremum cases eliminated)", titlefontsize=10)
    plot!(xlabel="true success probability", guidefontsize=9)
    plot!(xtick=0:0.1:1)
    plot!(legend=legend_KL, legendfontsize=7)
    plot!(true_prob, EE_KL_MLE.(true_prob);    label="MLE")
    plot!(true_prob, EE_KL_BayesJ.(true_prob); label="Bayes($β₁; $a₁, $b₁)", ls=:dash)
    plot!(true_prob, EE_KL_Bayes.(true_prob);  label="Bayes($β₂; $a₂, $b₂)", ls=:dashdot)
    
    P2 = plot(title="mean squared error (n = $n)", titlefontsize=10)
    plot!(xlabel="true success probability", guidefontsize=9)
    plot!(xtick=0:0.1:1)
    plot!(legend=legend_MSE, legendfontsize=7)
    plot!(true_prob, E_MSE_MLE.(true_prob);    label="MLE")
    plot!(true_prob, E_MSE_BayesJ.(true_prob); label="Bayes($β₁; $a₁, $b₁)", ls=:dash)
    plot!(true_prob, E_MSE_Bayes.(true_prob);  label="Bayes($β₂; $a₂, $b₂)", ls=:dashdot)
    
    plot(P1, P2; size=(800, 250), layout=(1,2))
end

Out[16]:

plot_mean_estimation_errors (generic function with 1 method)

In [17]:

plot_mean_estimation_errors(10)

Out[17]:

In [18]:

plot_mean_estimation_errors(20)

Out[18]:

In [19]:

plot_mean_estimation_errors(50; legend_KL=:top)

Out[19]:

In [20]:

plot_mean_estimation_errors(100)

Out[20]:

AICとWAICとLOOCV¶

以下では, 渡辺澄夫著『ベイズ統計の理論と方法』におけるAICとWAICとLOOCV(1個抜き出し交差検証)の定義を aic, waic, loocv と書くことにする. それらの定義は伝統的なAICの定義の $2n$ 分の1になっている. 伝統的なスケールの側での AIC, WAIC, LOOCV を AIC, WAIC, LOOCV と書くことにする.

KL情報量で定義された予測誤差を pe (prediction error)と書き, その $2n$ 倍を PE と書くことにする. 予測分布の汎化誤差を ge と書き, その $2n$ 倍を GE と書くことにする.

In [21]:

# 予測分布の対数尤度の -1/n 倍
function T(n, k; β=1.0, a=0.5, b=0.5, w̄=(β*k+a)/(β*n+a+b))
    p = k/n
    -(safexlogy(p, w̄) + safexlogy(1-p, 1-w̄))
end

# 事後分布に関する log p(X_i|w) の分散のiに関する和の 1/n 倍
function V(n, k; β=1.0, a=0.5, b=0.5)
    p = k/n
    p*trigamma(β*k+a) + (1-p)*trigamma(β*(n-k)+b) - trigamma(β*n+a+b)
end

aic(n, k) = T(n, k; a=0.0, b=0.0) + 1/n
waic(n, k; β=1.0, a=0.5, b=0.5) = T(n, k; β=β, a=a, b=b) + β*V(n, k; β=β, a=a, b=b)
function loocv_old(n, k; β=1.0, a=0.5, b=0.5)
    p = k/n
    (
        -   p  *(logbeta((1-β)+β*k+a, β*(n-k)+b) - logbeta(-β+β*k+a, β*(n-k)+b))
        - (1-p)*(logbeta(β*k+a, (1-β)+β*(n-k)+b) - logbeta(β*k+a, -β+β*(n-k)+b))
    )
end
function loocv(n, k; β=1.0, a=0.5, b=0.5)
    p = k/n
    -safexlogy(p, -β+β*k+a) - safexlogy(1-p, -β+β*(n-k)+b) + log(-β+β*n+a+b)
end

AIC(n, k)  = 2n*aic(n, k)
WAIC(n, k; β=1.0, a=0.5, b=0.5) = 2n*waic(n, k; β=β, a=a, b=b)
LOOCV(n, k; β=1.0, a=0.5, b=0.5) = 2n*loocv(n, k; β=β, a=a, b=b)

# https://ja.wikipedia.org/wiki/%E8%B5%A4%E6%B1%A0%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96
# 
#    AICc = AIC + 2nparam(nparam + 1)/(samplesize - nparam - 1)
#
# nparam = 1,  samplesize = n
#
AICc(n, k) = AIC(n, k) + (2*1^2 + 2*1)/(n-1-1)

ge(pred_dist, q, n, k) = G(q, pred_dist(n, k))
pe(pred_dist, q, n, k) = KL(q, pred_dist(n, k))

GE(pred_dist, q, n, k) = 2n*ge(pred_dist, q, n, k)
PE(pred_dist, q, n, k) = 2n*pe(pred_dist, q, n, k)

aic0(q, n, k) = T(n, k; a=0.0, b=0.0, w̄=q.p)
waic0(q, n, k; β=1.0, a=0.5, b=0.5) = T(n, k; β=β, a=a, b=b, w̄=q.p)
loocv0(q, n, k; β=1.0, a=0.5, b=0.5) = waic0(q, n, k; β=β, a=a, b=b)

AIC0(q, n, k) = 2n*aic0(q, n, k)
WAIC0(q, n, k; β=1.0, a=0.5, b=0.5) = 2n*waic0(q, n, k; β=β, a=a, b=b)
LOOCV0(q, n, k; β=1.0, a=0.5, b=0.5) = 2n*loocv0(q, n, k; β=β, a=a, b=b)

function supp(dist, m=2)
    μ, σ = mean(dist), std(dist)
    a = max(minimum(dist), floor(Int, μ - m*σ))
    b = min(maximum(dist),  ceil(Int, μ + m*σ))
    a:b
end

sample_supp(q, n, m=2) = supp(Binomial(n, q.p), m)

E(q, n, f) = E(Binomial(n, q.p), f)
EE(q, n, f) = EE(Binomial(n, q.p), f)

Out[21]:

EE (generic function with 2 methods)

In [22]:

loocv_old(10, 0) # error

DomainError with -0.5:
`gamma(x)` must be non-negative

Stacktrace:
 [1] loggamma at C:\Users\genkuroki\.julia\packages\SpecialFunctions\O3MVw\src\gamma.jl:651 [inlined]
 [2] logbeta(::Float64, ::Float64) at C:\Users\genkuroki\.julia\packages\SpecialFunctions\O3MVw\src\gamma.jl:797
 [3] #loocv_old#69(::Float64, ::Float64, ::Float64, ::typeof(loocv_old), ::Int64, ::Int64) at .\In[21]:17
 [4] loocv_old(::Int64, ::Int64) at .\In[21]:16
 [5] top-level scope at In[22]:1

In [23]:

loocv(10, 0)

Out[23]:

0.05129329438755059

平均AICと平均WAICと平均LOOCVと平均汎化誤差の関係¶

$n$ を大きくすると、平均AICと平均WAICと平均LOOCVと平均汎化誤差はほぼ等しくなる. (ただし, ここで平均とはサンプルの標本分布に関する平均のことである.)

In [24]:

function plot_EEs_and_Es(n; f=[fill(true, 5); false])
    w = range(0.005, 0.995, step=0.005)
    q = Bernoulli.(w)
    EE_GE_MLE  = (q -> EE(q, n, k->GE(pred_MLE, q, n, k))).(q)
    EE_AIC     = (q -> EE(q, n, k->AIC(n, k))).(q)
    EE_AICc    = (q -> EE(q, n, k->AICc(n, k))).(q)
    E_GE_Bayes = (q -> E(q, n, k->GE(pred_Bayes, q, n, k))).(q)
    E_WAIC     = (q -> E(q, n, k->WAIC(n, k))).(q)
    E_LOOCV    = (q -> E(q, n, k->LOOCV(n, k))).(q)

    pal = palette(:default)
    P = plot(title="n = $n", titlefontsize=10)
    plot!(xlabel="w", ylabel="expectation value", xtick=0:0.1:1)
    f[1] && plot!(w, EE_GE_MLE;  label="EE[GE_MLE]",   ls=:dot,     color=pal[1])
    f[2] && plot!(w, EE_AIC;     label="EE[AIC]",      ls=:dash,    color=pal[2])
    f[3] && plot!(w, E_GE_Bayes; label="E[GE_Bayes]",  ls=:solid,   color=pal[3])
    f[4] && plot!(w, E_WAIC;     label="E[WAIC]",      ls=:dash,    color=pal[4])
    f[5] && plot!(w, E_LOOCV;    label="E[LOOCV]",     ls=:dashdot, color=pal[5])
    f[6] && plot!(w, EE_AICc;     label="EE[AICc]",    ls=:dash,    color=pal[6])
    
    plot(P; size=(400, 250))
end

Out[24]:

plot_EEs_and_Es (generic function with 1 method)

In [25]:

plot_EEs_and_Es(10)

Out[25]:

$n$ が小さいとき, 平均AICの誤差は大きい.

In [26]:

plot_EEs_and_Es(20)

Out[26]:

In [27]:

plot_EEs_and_Es(50)

Out[27]:

In [28]:

plot_EEs_and_Es(100)

Out[28]:

In [29]:

plot_EEs_and_Es(10; f=[true; true; fill(false, 3); true])

Out[29]:

In [30]:

plot_EEs_and_Es(20; f=[true; true; fill(false, 3); true])

Out[30]:

AICとWAICとLOOCVと予測誤差の関係¶

渡辺澄夫『ベイズ推定の理論と方法』のp.80の12行目の公式(の両辺を $2n$ 倍したもの)とp.119の定理15(の公式の両辺を $2n$ 倍したもの)の数値的確認. 我々はパラメーターの個数が $d=1$ の場合の正則モデルを扱っているので, $\lambda=d/2$ となり, $\beta = 1$ のとき, 右辺を $2n$ 倍した結果は $n$ を大きくすると, $2d=4\lambda=2$ に近付くはずである.

In [31]:

function make_xtick(n, w)
    q = Bernoulli(w)
    q_sample = Binomial(n, q.p)
    
    t3 = round(Int, mean(q_sample))
    t2, t4 = extrema(sample_supp(q, n, 1))
    t1, t5 = extrema(sample_supp(q, n, 2))
    xticklabel = fill("", t5 - t1 + 1)
    xticklabel[1] = "$t1"
    xticklabel[t2-t1+1] = "$t2"
    xticklabel[t3-t1+1] = "$t3"
    xticklabel[t4-t1+1] = "$t4"
    xticklabel[end] = "$t5"
    tmp = cdf.(q_sample, t1-1:t5)
    xtick = @views @. (tmp[1:end-1] + tmp[2:end])/2
    
    xtick, xticklabel
end

function plot_AIC(n, w;
        legend = :top,
        ylim = (-4, 6),
        ytick = ylim[1]:ylim[2]
    )

    q = Bernoulli(w)
    q_sample = Binomial(n, q.p)

    x = range(0, 1, length=801)
    k = quantile.(q_sample, x)
    xtick, xticklabel = make_xtick(n, w)
    
    pred_error = @. PE(pred_MLE, q, n, k)
    AIC_minus_AIC0 = @. AIC(n,k) - AIC0(q,n,k)
    
    P = plot(title="n = $n,  w = $w"; titlefontsize=10)
    plot!(xlim=(0, 1), ylim=ylim)
    plot!(xtick=(xtick, xticklabel), ytick=ytick)
    plot!(xlabel="k", guidefontsize=10)
    plot!(legend=legend)
    plot!(x, pred_error; label="Prediction Error of MLE")
    plot!(x, AIC_minus_AIC0; label="AIC - AIC0", ls=:dash, lw=1.3)
    hline!([0]; color=:black, ls=:dot, alpha=0.8, label="")
    
    plot(P; size=(500, 350))
end

function plot_WAIC_LOOCV(n, w;
        legend = :top,
        ylim = (-4, 6),
        ytick = ylim[1]:ylim[2]
    )

    q = Bernoulli(w)
    q_sample = Binomial(n, w)

    x = range(0, 1, length=801)
    k = quantile.(q_sample, x)
    xtick, xticklabel = make_xtick(n, w)
    
    pred_error = @. PE(pred_Bayes, q, n, k)
    WAIC_minus_WAIC0 = @. WAIC(n,k) - WAIC0(q,n,k)
    LOOCV_minus_LOOCV0 = @. LOOCV(n,k) - LOOCV0(q,n,k)
    
    P = plot(title="n = $n,  w = $w"; titlefontsize=10)
    plot!(xlim=(0, 1), ylim=ylim)
    plot!(xtick=(xtick, xticklabel), ytick=ytick)
    plot!(xlabel="k", guidefontsize=10)
    plot!(legend=legend)
    plot!(x, pred_error; label="Prediction Error of Bayes")
    plot!(x, WAIC_minus_WAIC0; label="WAIC - WAIC0", ls=:dash, lw=1.3)
    plot!(x, LOOCV_minus_LOOCV0; label="LOOCV - LOOCV0", ls=:dot, lw=1.3)
    hline!([0]; color=:black, ls=:dot, alpha=0.8, label="")
    
    plot(P; size=(500, 350))
end

Out[31]:

plot_WAIC_LOOCV (generic function with 1 method)

In [32]:

n = 10
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_AIC = E(q, n, k->AIC(n,k) < AIC0(q,n,k)) |> rd
@show prob_fail_AIC
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_MLE, q, n, k)
    AIC_minus_AIC0 = AIC(n,k) - AIC0(q,n,k)
    (
        k = k,
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_MLE = pred_error |>rd, 
        AIC_minus_AIC0 = AIC_minus_AIC0 |>rd, 
        PE_plus_AIC_minus_AIC0 = pred_error + AIC_minus_AIC0 |>rd
    ) |> println
end

n = 10
w = 0.4
prob_fail_AIC = 0.101

(k = 0, prob = 0.006, PE_MLE = Inf, AIC_minus_AIC0 = -8.217, PE_plus_AIC_minus_AIC0 = Inf)
(k = 1, prob = 0.04, PE_MLE = 6.225, AIC_minus_AIC0 = -2.526, PE_plus_AIC_minus_AIC0 = 3.699)
(k = 2, prob = 0.121, PE_MLE = 2.093, AIC_minus_AIC0 = 0.17, PE_plus_AIC_minus_AIC0 = 2.263)
(k = 3, prob = 0.215, PE_MLE = 0.452, AIC_minus_AIC0 = 1.568, PE_plus_AIC_minus_AIC0 = 2.02)
(k = 4, prob = 0.251, PE_MLE = 0.0, AIC_minus_AIC0 = 2.0, PE_plus_AIC_minus_AIC0 = 2.0)
(k = 5, prob = 0.201, PE_MLE = 0.403, AIC_minus_AIC0 = 1.592, PE_plus_AIC_minus_AIC0 = 1.994)
(k = 6, prob = 0.111, PE_MLE = 1.622, AIC_minus_AIC0 = 0.378, PE_plus_AIC_minus_AIC0 = 2.0)
(k = 7, prob = 0.042, PE_MLE = 3.841, AIC_minus_AIC0 = -1.676, PE_plus_AIC_minus_AIC0 = 2.165)
(k = 8, prob = 0.011, PE_MLE = 7.638, AIC_minus_AIC0 = -4.696, PE_plus_AIC_minus_AIC0 = 2.942)
(k = 9, prob = 0.002, PE_MLE = 15.014, AIC_minus_AIC0 = -9.013, PE_plus_AIC_minus_AIC0 = 6.0)

In [33]:

n = 10
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_WAIC = E(q, n, k -> WAIC(n,k) < WAIC0(q,n,k)) |> rd
@show prob_fail_WAIC
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_Bayes, q, n, k)
    WAIC_minus_WAIC0 = WAIC(n,k) - WAIC0(q,n,k)
    (
        k = k, 
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_Bayes = pred_error |>rd, 
        WAIC_minus_WAIC0 = WAIC_minus_WAIC0 |>rd, 
        PE_plus_WAIC_minus_WAIC0 = pred_error + WAIC_minus_WAIC0 |>rd
    ) |> println
end

n = 10
w = 0.4
prob_fail_WAIC = 0.101

(k = 0, prob = 0.006, PE_Bayes = 11.826, WAIC_minus_WAIC0 = -9.191, PE_plus_WAIC_minus_WAIC0 = 2.635)
(k = 1, prob = 0.04, PE_Bayes = 4.238, WAIC_minus_WAIC0 = -2.439, PE_plus_WAIC_minus_WAIC0 = 1.799)
(k = 2, prob = 0.121, PE_Bayes = 1.487, WAIC_minus_WAIC0 = 0.269, PE_plus_WAIC_minus_WAIC0 = 1.755)
(k = 3, prob = 0.215, PE_Bayes = 0.297, WAIC_minus_WAIC0 = 1.659, PE_plus_WAIC_minus_WAIC0 = 1.956)
(k = 4, prob = 0.251, PE_Bayes = 0.003, WAIC_minus_WAIC0 = 2.085, PE_plus_WAIC_minus_WAIC0 = 2.089)
(k = 5, prob = 0.201, PE_Bayes = 0.403, WAIC_minus_WAIC0 = 1.675, PE_plus_WAIC_minus_WAIC0 = 2.078)
(k = 6, prob = 0.111, PE_Bayes = 1.474, WAIC_minus_WAIC0 = 0.463, PE_plus_WAIC_minus_WAIC0 = 1.938)
(k = 7, prob = 0.042, PE_Bayes = 3.345, WAIC_minus_WAIC0 = -1.585, PE_plus_WAIC_minus_WAIC0 = 1.76)
(k = 8, prob = 0.011, PE_Bayes = 6.382, WAIC_minus_WAIC0 = -4.597, PE_plus_WAIC_minus_WAIC0 = 1.785)
(k = 9, prob = 0.002, PE_Bayes = 11.622, WAIC_minus_WAIC0 = -8.927, PE_plus_WAIC_minus_WAIC0 = 2.695)

In [34]:

n = 10
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_LOOCV = EE(q, n, k->LOOCV(n,k) < LOOCV0(q,n,k)) |> rd
@show prob_fail_LOOCV
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_Bayes, q, n, k)
    LOOCV_minus_LOOCV0 = LOOCV(n,k) - LOOCV0(q,n,k)
    (
        k = k, 
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_Bayes = pred_error |>rd, 
        LOOCV_minus_LOOCVC0 = LOOCV_minus_LOOCV0 |>rd, 
        PE_plus_LOOCV_minus_LOOCV0 = pred_error + LOOCV_minus_LOOCV0 |>rd
    ) |> println
end

n = 10
w = 0.4
prob_fail_LOOCV = 0.096

(k = 0, prob = 0.006, PE_Bayes = 11.826, LOOCV_minus_LOOCVC0 = -9.191, PE_plus_LOOCV_minus_LOOCV0 = 2.636)
(k = 1, prob = 0.04, PE_Bayes = 4.238, LOOCV_minus_LOOCVC0 = -2.111, PE_plus_LOOCV_minus_LOOCV0 = 2.128)
(k = 2, prob = 0.121, PE_Bayes = 1.487, LOOCV_minus_LOOCVC0 = 0.353, PE_plus_LOOCV_minus_LOOCV0 = 1.84)
(k = 3, prob = 0.215, PE_Bayes = 0.297, LOOCV_minus_LOOCVC0 = 1.699, PE_plus_LOOCV_minus_LOOCV0 = 1.996)
(k = 4, prob = 0.251, PE_Bayes = 0.003, LOOCV_minus_LOOCVC0 = 2.112, PE_plus_LOOCV_minus_LOOCV0 = 2.116)
(k = 5, prob = 0.201, PE_Bayes = 0.403, LOOCV_minus_LOOCVC0 = 1.699, PE_plus_LOOCV_minus_LOOCV0 = 2.102)
(k = 6, prob = 0.111, PE_Bayes = 1.474, LOOCV_minus_LOOCVC0 = 0.491, PE_plus_LOOCV_minus_LOOCV0 = 1.965)
(k = 7, prob = 0.042, PE_Bayes = 3.345, LOOCV_minus_LOOCVC0 = -1.544, PE_plus_LOOCV_minus_LOOCV0 = 1.801)
(k = 8, prob = 0.011, PE_Bayes = 6.382, LOOCV_minus_LOOCVC0 = -4.513, PE_plus_LOOCV_minus_LOOCV0 = 1.869)
(k = 9, prob = 0.002, PE_Bayes = 11.622, LOOCV_minus_LOOCVC0 = -8.598, PE_plus_LOOCV_minus_LOOCV0 = 3.024)

In [35]:

n = 10
w = 0.4
plot_AIC(n, w) |> display
plot_WAIC_LOOCV(n, w) |> display
pyplotclf()

In [36]:

n = 20
w = 0.4
plot_AIC(n, w) |> display
plot_WAIC_LOOCV(n, w) |> display
pyplotclf()

In [37]:

n = 100
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_AIC = E(q, n, k->AIC(n,k) < AIC0(q,n,k)) |> rd
@show prob_fail_AIC
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_MLE, q, n, k)
    AIC_minus_AIC0 = AIC(n,k) - AIC0(q,n,k)
    (
        k = k,
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_MLE = pred_error |>rd, 
        AIC_minus_AIC0 = AIC_minus_AIC0 |>rd, 
        PE_plus_AIC_minus_AIC0 = pred_error + AIC_minus_AIC0 |>rd
    ) |> println
end

n = 100
w = 0.4
prob_fail_AIC = 0.184

(k = 25, prob = 0.001, PE_MLE = 10.823, AIC_minus_AIC0 = -7.971, PE_plus_AIC_minus_AIC0 = 2.852)
(k = 26, prob = 0.001, PE_MLE = 9.296, AIC_minus_AIC0 = -6.638, PE_plus_AIC_minus_AIC0 = 2.658)
(k = 27, prob = 0.002, PE_MLE = 7.91, AIC_minus_AIC0 = -5.408, PE_plus_AIC_minus_AIC0 = 2.501)
(k = 28, prob = 0.004, PE_MLE = 6.655, AIC_minus_AIC0 = -4.281, PE_plus_AIC_minus_AIC0 = 2.375)
(k = 29, prob = 0.006, PE_MLE = 5.526, AIC_minus_AIC0 = -3.252, PE_plus_AIC_minus_AIC0 = 2.275)
(k = 30, prob = 0.01, PE_MLE = 4.516, AIC_minus_AIC0 = -2.32, PE_plus_AIC_minus_AIC0 = 2.196)
(k = 31, prob = 0.015, PE_MLE = 3.62, AIC_minus_AIC0 = -1.484, PE_plus_AIC_minus_AIC0 = 2.136)
(k = 32, prob = 0.022, PE_MLE = 2.832, AIC_minus_AIC0 = -0.741, PE_plus_AIC_minus_AIC0 = 2.091)
(k = 33, prob = 0.03, PE_MLE = 2.148, AIC_minus_AIC0 = -0.09, PE_plus_AIC_minus_AIC0 = 2.058)
(k = 34, prob = 0.039, PE_MLE = 1.564, AIC_minus_AIC0 = 0.47, PE_plus_AIC_minus_AIC0 = 2.035)
(k = 35, prob = 0.049, PE_MLE = 1.077, AIC_minus_AIC0 = 0.942, PE_plus_AIC_minus_AIC0 = 2.019)
(k = 36, prob = 0.059, PE_MLE = 0.684, AIC_minus_AIC0 = 1.325, PE_plus_AIC_minus_AIC0 = 2.009)
(k = 37, prob = 0.068, PE_MLE = 0.382, AIC_minus_AIC0 = 1.622, PE_plus_AIC_minus_AIC0 = 2.004)
(k = 38, prob = 0.075, PE_MLE = 0.169, AIC_minus_AIC0 = 1.832, PE_plus_AIC_minus_AIC0 = 2.001)
(k = 39, prob = 0.08, PE_MLE = 0.042, AIC_minus_AIC0 = 1.958, PE_plus_AIC_minus_AIC0 = 2.0)
(k = 40, prob = 0.081, PE_MLE = 0.0, AIC_minus_AIC0 = 2.0, PE_plus_AIC_minus_AIC0 = 2.0)
(k = 41, prob = 0.079, PE_MLE = 0.041, AIC_minus_AIC0 = 1.958, PE_plus_AIC_minus_AIC0 = 2.0)
(k = 42, prob = 0.074, PE_MLE = 0.165, AIC_minus_AIC0 = 1.834, PE_plus_AIC_minus_AIC0 = 1.999)
(k = 43, prob = 0.067, PE_MLE = 0.37, AIC_minus_AIC0 = 1.628, PE_plus_AIC_minus_AIC0 = 1.997)
(k = 44, prob = 0.058, PE_MLE = 0.654, AIC_minus_AIC0 = 1.34, PE_plus_AIC_minus_AIC0 = 1.994)
(k = 45, prob = 0.048, PE_MLE = 1.019, AIC_minus_AIC0 = 0.971, PE_plus_AIC_minus_AIC0 = 1.99)
(k = 46, prob = 0.038, PE_MLE = 1.462, AIC_minus_AIC0 = 0.521, PE_plus_AIC_minus_AIC0 = 1.983)
(k = 47, prob = 0.029, PE_MLE = 1.985, AIC_minus_AIC0 = -0.01, PE_plus_AIC_minus_AIC0 = 1.975)
(k = 48, prob = 0.021, PE_MLE = 2.586, AIC_minus_AIC0 = -0.62, PE_plus_AIC_minus_AIC0 = 1.966)
(k = 49, prob = 0.015, PE_MLE = 3.267, AIC_minus_AIC0 = -1.311, PE_plus_AIC_minus_AIC0 = 1.956)
(k = 50, prob = 0.01, PE_MLE = 4.027, AIC_minus_AIC0 = -2.082, PE_plus_AIC_minus_AIC0 = 1.945)
(k = 51, prob = 0.007, PE_MLE = 4.867, AIC_minus_AIC0 = -2.933, PE_plus_AIC_minus_AIC0 = 1.934)
(k = 52, prob = 0.004, PE_MLE = 5.788, AIC_minus_AIC0 = -3.864, PE_plus_AIC_minus_AIC0 = 1.924)
(k = 53, prob = 0.003, PE_MLE = 6.791, AIC_minus_AIC0 = -4.875, PE_plus_AIC_minus_AIC0 = 1.915)
(k = 54, prob = 0.001, PE_MLE = 7.876, AIC_minus_AIC0 = -5.967, PE_plus_AIC_minus_AIC0 = 1.909)
(k = 55, prob = 0.001, PE_MLE = 9.046, AIC_minus_AIC0 = -7.139, PE_plus_AIC_minus_AIC0 = 1.907)

In [38]:

n = 100
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_WAIC = E(q, n, k -> WAIC(n,k) < WAIC0(q,n,k)) |> rd
@show prob_fail_WAIC
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_Bayes, q, n, k)
    WAIC_minus_WAIC0 = WAIC(n,k) - WAIC0(q,n,k)
    (
        k = k, 
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_Bayes = pred_error |>rd, 
        WAIC_minus_WAIC0 = WAIC_minus_WAIC0 |>rd, 
        PE_plus_WAIC_minus_WAIC0 = pred_error + WAIC_minus_WAIC0 |>rd
    ) |> println
end

n = 100
w = 0.4
prob_fail_WAIC = 0.155

(k = 25, prob = 0.001, PE_Bayes = 10.432, WAIC_minus_WAIC0 = -7.958, PE_plus_WAIC_minus_WAIC0 = 2.473)
(k = 26, prob = 0.001, PE_Bayes = 8.954, WAIC_minus_WAIC0 = -6.625, PE_plus_WAIC_minus_WAIC0 = 2.329)
(k = 27, prob = 0.002, PE_Bayes = 7.613, WAIC_minus_WAIC0 = -5.396, PE_plus_WAIC_minus_WAIC0 = 2.216)
(k = 28, prob = 0.004, PE_Bayes = 6.399, WAIC_minus_WAIC0 = -4.268, PE_plus_WAIC_minus_WAIC0 = 2.131)
(k = 29, prob = 0.006, PE_Bayes = 5.307, WAIC_minus_WAIC0 = -3.24, PE_plus_WAIC_minus_WAIC0 = 2.067)
(k = 30, prob = 0.01, PE_Bayes = 4.33, WAIC_minus_WAIC0 = -2.309, PE_plus_WAIC_minus_WAIC0 = 2.022)
(k = 31, prob = 0.015, PE_Bayes = 3.464, WAIC_minus_WAIC0 = -1.472, PE_plus_WAIC_minus_WAIC0 = 1.991)
(k = 32, prob = 0.022, PE_Bayes = 2.703, WAIC_minus_WAIC0 = -0.73, PE_plus_WAIC_minus_WAIC0 = 1.973)
(k = 33, prob = 0.03, PE_Bayes = 2.043, WAIC_minus_WAIC0 = -0.079, PE_plus_WAIC_minus_WAIC0 = 1.964)
(k = 34, prob = 0.039, PE_Bayes = 1.481, WAIC_minus_WAIC0 = 0.481, PE_plus_WAIC_minus_WAIC0 = 1.962)
(k = 35, prob = 0.049, PE_Bayes = 1.013, WAIC_minus_WAIC0 = 0.952, PE_plus_WAIC_minus_WAIC0 = 1.966)
(k = 36, prob = 0.059, PE_Bayes = 0.637, WAIC_minus_WAIC0 = 1.336, PE_plus_WAIC_minus_WAIC0 = 1.973)
(k = 37, prob = 0.068, PE_Bayes = 0.35, WAIC_minus_WAIC0 = 1.632, PE_plus_WAIC_minus_WAIC0 = 1.982)
(k = 38, prob = 0.075, PE_Bayes = 0.149, WAIC_minus_WAIC0 = 1.843, PE_plus_WAIC_minus_WAIC0 = 1.992)
(k = 39, prob = 0.08, PE_Bayes = 0.033, WAIC_minus_WAIC0 = 1.969, PE_plus_WAIC_minus_WAIC0 = 2.002)
(k = 40, prob = 0.081, PE_Bayes = 0.0, WAIC_minus_WAIC0 = 2.01, PE_plus_WAIC_minus_WAIC0 = 2.011)
(k = 41, prob = 0.079, PE_Bayes = 0.049, WAIC_minus_WAIC0 = 1.969, PE_plus_WAIC_minus_WAIC0 = 2.018)
(k = 42, prob = 0.074, PE_Bayes = 0.178, WAIC_minus_WAIC0 = 1.844, PE_plus_WAIC_minus_WAIC0 = 2.023)
(k = 43, prob = 0.067, PE_Bayes = 0.387, WAIC_minus_WAIC0 = 1.638, PE_plus_WAIC_minus_WAIC0 = 2.025)
(k = 44, prob = 0.058, PE_Bayes = 0.674, WAIC_minus_WAIC0 = 1.35, PE_plus_WAIC_minus_WAIC0 = 2.024)
(k = 45, prob = 0.048, PE_Bayes = 1.039, WAIC_minus_WAIC0 = 0.981, PE_plus_WAIC_minus_WAIC0 = 2.02)
(k = 46, prob = 0.038, PE_Bayes = 1.482, WAIC_minus_WAIC0 = 0.531, PE_plus_WAIC_minus_WAIC0 = 2.012)
(k = 47, prob = 0.029, PE_Bayes = 2.002, WAIC_minus_WAIC0 = 0.0, PE_plus_WAIC_minus_WAIC0 = 2.002)
(k = 48, prob = 0.021, PE_Bayes = 2.599, WAIC_minus_WAIC0 = -0.611, PE_plus_WAIC_minus_WAIC0 = 1.989)
(k = 49, prob = 0.015, PE_Bayes = 3.274, WAIC_minus_WAIC0 = -1.301, PE_plus_WAIC_minus_WAIC0 = 1.973)
(k = 50, prob = 0.01, PE_Bayes = 4.027, WAIC_minus_WAIC0 = -2.072, PE_plus_WAIC_minus_WAIC0 = 1.955)
(k = 51, prob = 0.007, PE_Bayes = 4.859, WAIC_minus_WAIC0 = -2.923, PE_plus_WAIC_minus_WAIC0 = 1.935)
(k = 52, prob = 0.004, PE_Bayes = 5.769, WAIC_minus_WAIC0 = -3.854, PE_plus_WAIC_minus_WAIC0 = 1.915)
(k = 53, prob = 0.003, PE_Bayes = 6.76, WAIC_minus_WAIC0 = -4.865, PE_plus_WAIC_minus_WAIC0 = 1.894)
(k = 54, prob = 0.001, PE_Bayes = 7.831, WAIC_minus_WAIC0 = -5.957, PE_plus_WAIC_minus_WAIC0 = 1.875)
(k = 55, prob = 0.001, PE_Bayes = 8.986, WAIC_minus_WAIC0 = -7.129, PE_plus_WAIC_minus_WAIC0 = 1.857)

In [39]:

n = 100
w = 0.4
@show n
@show w
q = Bernoulli(w)
prob_fail_LOOCV = EE(q, n, k->LOOCV(n,k) < LOOCV0(q,n,k)) |> rd
@show prob_fail_LOOCV
println()
for k in sample_supp(q, n, 3)
    pred_error = PE(pred_Bayes, q, n, k)
    LOOCV_minus_LOOCV0 = LOOCV(n,k) - LOOCV0(q,n,k)
    (
        k = k, 
        prob = pdf(Binomial(n, q.p), k)|>rd, 
        PE_Bayes = pred_error |>rd, 
        LOOCV_minus_LOOCVC0 = LOOCV_minus_LOOCV0 |>rd, 
        PE_plus_LOOCV_minus_LOOCV0 = pred_error + LOOCV_minus_LOOCV0 |>rd
    ) |> println
end

n = 100
w = 0.4
prob_fail_LOOCV = 0.155

(k = 25, prob = 0.001, PE_Bayes = 10.432, LOOCV_minus_LOOCVC0 = -7.958, PE_plus_LOOCV_minus_LOOCV0 = 2.474)
(k = 26, prob = 0.001, PE_Bayes = 8.954, LOOCV_minus_LOOCVC0 = -6.625, PE_plus_LOOCV_minus_LOOCV0 = 2.33)
(k = 27, prob = 0.002, PE_Bayes = 7.613, LOOCV_minus_LOOCVC0 = -5.396, PE_plus_LOOCV_minus_LOOCV0 = 2.217)
(k = 28, prob = 0.004, PE_Bayes = 6.399, LOOCV_minus_LOOCVC0 = -4.268, PE_plus_LOOCV_minus_LOOCV0 = 2.131)
(k = 29, prob = 0.006, PE_Bayes = 5.307, LOOCV_minus_LOOCVC0 = -3.24, PE_plus_LOOCV_minus_LOOCV0 = 2.067)
(k = 30, prob = 0.01, PE_Bayes = 4.33, LOOCV_minus_LOOCVC0 = -2.308, PE_plus_LOOCV_minus_LOOCV0 = 2.022)
(k = 31, prob = 0.015, PE_Bayes = 3.464, LOOCV_minus_LOOCVC0 = -1.472, PE_plus_LOOCV_minus_LOOCV0 = 1.992)
(k = 32, prob = 0.022, PE_Bayes = 2.703, LOOCV_minus_LOOCVC0 = -0.729, PE_plus_LOOCV_minus_LOOCV0 = 1.973)
(k = 33, prob = 0.03, PE_Bayes = 2.043, LOOCV_minus_LOOCVC0 = -0.079, PE_plus_LOOCV_minus_LOOCV0 = 1.964)
(k = 34, prob = 0.039, PE_Bayes = 1.481, LOOCV_minus_LOOCVC0 = 0.482, PE_plus_LOOCV_minus_LOOCV0 = 1.962)
(k = 35, prob = 0.049, PE_Bayes = 1.013, LOOCV_minus_LOOCVC0 = 0.953, PE_plus_LOOCV_minus_LOOCV0 = 1.966)
(k = 36, prob = 0.059, PE_Bayes = 0.637, LOOCV_minus_LOOCVC0 = 1.336, PE_plus_LOOCV_minus_LOOCV0 = 1.973)
(k = 37, prob = 0.068, PE_Bayes = 0.35, LOOCV_minus_LOOCVC0 = 1.632, PE_plus_LOOCV_minus_LOOCV0 = 1.982)
(k = 38, prob = 0.075, PE_Bayes = 0.149, LOOCV_minus_LOOCVC0 = 1.843, PE_plus_LOOCV_minus_LOOCV0 = 1.992)
(k = 39, prob = 0.08, PE_Bayes = 0.033, LOOCV_minus_LOOCVC0 = 1.969, PE_plus_LOOCV_minus_LOOCV0 = 2.002)
(k = 40, prob = 0.081, PE_Bayes = 0.0, LOOCV_minus_LOOCVC0 = 2.01, PE_plus_LOOCV_minus_LOOCV0 = 2.011)
(k = 41, prob = 0.079, PE_Bayes = 0.049, LOOCV_minus_LOOCVC0 = 1.969, PE_plus_LOOCV_minus_LOOCV0 = 2.018)
(k = 42, prob = 0.074, PE_Bayes = 0.178, LOOCV_minus_LOOCVC0 = 1.845, PE_plus_LOOCV_minus_LOOCV0 = 2.023)
(k = 43, prob = 0.067, PE_Bayes = 0.387, LOOCV_minus_LOOCVC0 = 1.638, PE_plus_LOOCV_minus_LOOCV0 = 2.025)
(k = 44, prob = 0.058, PE_Bayes = 0.674, LOOCV_minus_LOOCVC0 = 1.35, PE_plus_LOOCV_minus_LOOCV0 = 2.024)
(k = 45, prob = 0.048, PE_Bayes = 1.039, LOOCV_minus_LOOCVC0 = 0.981, PE_plus_LOOCV_minus_LOOCV0 = 2.02)
(k = 46, prob = 0.038, PE_Bayes = 1.482, LOOCV_minus_LOOCVC0 = 0.531, PE_plus_LOOCV_minus_LOOCV0 = 2.012)
(k = 47, prob = 0.029, PE_Bayes = 2.002, LOOCV_minus_LOOCVC0 = 0.0, PE_plus_LOOCV_minus_LOOCV0 = 2.002)
(k = 48, prob = 0.021, PE_Bayes = 2.599, LOOCV_minus_LOOCVC0 = -0.61, PE_plus_LOOCV_minus_LOOCV0 = 1.989)
(k = 49, prob = 0.015, PE_Bayes = 3.274, LOOCV_minus_LOOCVC0 = -1.301, PE_plus_LOOCV_minus_LOOCV0 = 1.973)
(k = 50, prob = 0.01, PE_Bayes = 4.027, LOOCV_minus_LOOCVC0 = -2.072, PE_plus_LOOCV_minus_LOOCV0 = 1.955)
(k = 51, prob = 0.007, PE_Bayes = 4.859, LOOCV_minus_LOOCVC0 = -2.923, PE_plus_LOOCV_minus_LOOCV0 = 1.935)
(k = 52, prob = 0.004, PE_Bayes = 5.769, LOOCV_minus_LOOCVC0 = -3.854, PE_plus_LOOCV_minus_LOOCV0 = 1.915)
(k = 53, prob = 0.003, PE_Bayes = 6.76, LOOCV_minus_LOOCVC0 = -4.865, PE_plus_LOOCV_minus_LOOCV0 = 1.895)
(k = 54, prob = 0.001, PE_Bayes = 7.831, LOOCV_minus_LOOCVC0 = -5.956, PE_plus_LOOCV_minus_LOOCV0 = 1.875)
(k = 55, prob = 0.001, PE_Bayes = 8.986, LOOCV_minus_LOOCVC0 = -7.128, PE_plus_LOOCV_minus_LOOCV0 = 1.857)

In [40]:

n = 100
w = 0.4
plot_AIC(n, w) |> display
plot_WAIC_LOOCV(n, w) |> display
pyplotclf()

In [41]:

n = 1000
w = 0.4
plot_AIC(n, w) |> display
plot_WAIC_LOOCV(n, w) |> display
pyplotclf()

AICとWAICとLOOCVによって正しいモデル選択に失敗する確率¶

In [42]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_AIC = E(q, n, k -> AIC(n,k) < AIC0(q,n,k))
    (n=n, prob_fail_AIC=prob_fail_AIC|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_AIC = 0.101119)
(n = 32, prob_fail_AIC = 0.149451)
(n = 100, prob_fail_AIC = 0.184234)
(n = 316, prob_fail_AIC = 0.168147)
(n = 1000, prob_fail_AIC = 0.165153)
(n = 3162, prob_fail_AIC = 0.156863)
(n = 10000, prob_fail_AIC = 0.155993)
(n = 31623, prob_fail_AIC = 0.157987)
(n = 100000, prob_fail_AIC = 0.156522)
(n = 316228, prob_fail_AIC = 0.157407)
(n = 1000000, prob_fail_AIC = 0.157491)

In [43]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_WAIC = E(q, n, k -> WAIC(n,k) < WAIC0(q,n,k))
    (n=n, prob_fail_WAIC=prob_fail_WAIC|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_WAIC = 0.101119)
(n = 32, prob_fail_WAIC = 0.149451)
(n = 100, prob_fail_WAIC = 0.155043)
(n = 316, prob_fail_WAIC = 0.168147)
(n = 1000, prob_fail_WAIC = 0.165153)
(n = 3162, prob_fail_WAIC = 0.156863)
(n = 10000, prob_fail_WAIC = 0.155993)
(n = 31623, prob_fail_WAIC = 0.157987)
(n = 100000, prob_fail_WAIC = 0.156522)
(n = 316228, prob_fail_WAIC = 0.157407)
(n = 1000000, prob_fail_WAIC = 0.157491)

In [44]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_LOOCV = EE(q, n, k -> LOOCV(n,k) < LOOCV0(q,n,k))
    (n=n, prob_fail_LOOCV=prob_fail_LOOCV|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_LOOCV = 0.095556)
(n = 32, prob_fail_LOOCV = 0.149451)
(n = 100, prob_fail_LOOCV = 0.155043)
(n = 316, prob_fail_LOOCV = 0.168147)
(n = 1000, prob_fail_LOOCV = 0.165153)
(n = 3162, prob_fail_LOOCV = 0.156863)
(n = 10000, prob_fail_LOOCV = 0.155993)
(n = 31623, prob_fail_LOOCV = 0.157987)
(n = 100000, prob_fail_LOOCV = 0.156522)
(n = 316228, prob_fail_LOOCV = 0.157407)
(n = 1000000, prob_fail_LOOCV = 0.157491)

In [45]:

ccdf(Chisq(1), 2)|>x->rd(x,6)

Out[45]:

0.157299

この場合には正しいモデル選択に失敗する確率は, 自由度1のχ²分布で2より大きくなる確率 0.157 に収束する. 対数尤度比検定との関係でこれは当然こうなるべきである.

BICとWBICと自由エネルギー¶

In [46]:

bic(n, k) = T(n, k; a=0.0, b=0.0) + (1/2)*log(n)/n
bic0(q, n, k) = aic0(q, n, k)

BIC(n, k) = 2n*bic(n, k)
BIC0(q, n, k) = 2n*bic0(q, n, k)

function wbic(n, k; a=0.5, b=0.5)
    β = 1/log(n)
    p = k/n
    -p*digamma(β*k+a) - (1-p)*digamma(β*(n-k)+b) + digamma(β*n+a+b)
end
wbic0(q, n, k; a=0.5, b=0.5) = waic0(q, n, k; β=1.0, a=a, b=b)

WBIC(n, k; a=0.5, b=0.5) = 2n*wbic(n, k; a=a, b=b)
WBIC0(q, n, k; a=0.5, b=0.5) = 2n*wbic0(q, n, k; a=a, b=b)

FE(n, k; β=1.0, a=0.5, b=0.5) = (2/β)*(
    - (loggamma(a+k) - loggamma(a)) 
    - (loggamma(b+n-k) - loggamma(b))
    + (loggamma(a+b+n) - loggamma(a+b))
)
FE0(q, n, k; β=1.0) = (2/β)*(- safemul(k, log(q.p)) - safemul(n-k, log(1-q.p)))

Out[46]:

FE0 (generic function with 1 method)

BICとWBICと自由エネルギーの関係¶

In [47]:

function plot_BIC_WBIC_FE(n, w;
        legend = :bottom,
        ylim = (-4, 6),
        ytick = ylim[1]:ylim[2]
    )

    q = Bernoulli(w)
    q_sample = Binomial(n, q.p)

    x = range(0, 1, length=801)
    k = quantile.(q_sample, x)
    xtick, xticklabel = make_xtick(n, w)
    
    FE_minus_FE0 = @. FE(n,k) - FE0(q,n,k)
    BIC_minus_BIC0 = @. BIC(n,k) - BIC0(q,n,k)
    WBIC_minus_WBIC0 = @. WBIC(n,k) - WBIC0(q,n,k)
    
    P = plot(title="n = $n,  w = $w"; titlefontsize=10)
    plot!(xlim=(0, 1), ylim=ylim)
    plot!(xtick=(xtick, xticklabel), ytick=ytick)
    plot!(xlabel="k", guidefontsize=10)
    plot!(legend=legend)
    plot!(x, FE_minus_FE0; label="FE - FE0")
    plot!(x, WBIC_minus_WBIC0; label="WBIC - WBIC0", ls=:dash, lw=1.2)
    plot!(x, BIC_minus_BIC0; label="BIC - BIC0", ls=:dot, lw=1.5)
    hline!([0]; color=:black, ls=:dot, alpha=0.8, label="")
    
    plot(P; size=(500, 350))
end

Out[47]:

plot_BIC_WBIC_FE (generic function with 1 method)

In [48]:

n = 10
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>rd, 
        FE0=FE0(q,n,k)|>rd, 
        BIC0=BIC0(q,n,k)|>rd, 
        WBIC0=WBIC0(q,n,k)|>rd, 
    ) |> println
end

n = 10
w = 0.4
(k = 0, prob = 0.006, FE0 = 10.217, BIC0 = 10.217, WBIC0 = 10.217)
(k = 1, prob = 0.04, FE0 = 11.027, BIC0 = 11.027, WBIC0 = 11.027)
(k = 2, prob = 0.121, FE0 = 11.838, BIC0 = 11.838, WBIC0 = 11.838)
(k = 3, prob = 0.215, FE0 = 12.649, BIC0 = 12.649, WBIC0 = 12.649)
(k = 4, prob = 0.251, FE0 = 13.46, BIC0 = 13.46, WBIC0 = 13.46)
(k = 5, prob = 0.201, FE0 = 14.271, BIC0 = 14.271, WBIC0 = 14.271)
(k = 6, prob = 0.111, FE0 = 15.082, BIC0 = 15.082, WBIC0 = 15.082)
(k = 7, prob = 0.042, FE0 = 15.893, BIC0 = 15.893, WBIC0 = 15.893)
(k = 8, prob = 0.011, FE0 = 16.704, BIC0 = 16.704, WBIC0 = 16.704)
(k = 9, prob = 0.002, FE0 = 17.515, BIC0 = 17.515, WBIC0 = 17.515)

In [49]:

n = 10
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    F0 = FE0(q,n,k)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>x->rd(x,4), 
        FE_minus_FE0=FE(n,k)-F0|>rd, 
        BIC_minus_FE0=BIC(n,k)-F0|>rd,
        WBIC_minus_FE0=WBIC(n,k)-F0|>rd
    ) |> println
end

plot_BIC_WBIC_FE(n, w)

n = 10
w = 0.4
(k = 0, prob = 0.006, FE_minus_FE0 = -6.744, BIC_minus_FE0 = -7.914, WBIC_minus_FE0 = -8.046)
(k = 1, prob = 0.0403, FE_minus_FE0 = -1.666, BIC_minus_FE0 = -2.223, WBIC_minus_FE0 = -2.646)
(k = 2, prob = 0.1209, FE_minus_FE0 = 0.992, BIC_minus_FE0 = 0.472, WBIC_minus_FE0 = 0.143)
(k = 3, prob = 0.215, FE_minus_FE0 = 2.378, BIC_minus_FE0 = 1.871, WBIC_minus_FE0 = 1.585)
(k = 4, prob = 0.2508, FE_minus_FE0 = 2.805, BIC_minus_FE0 = 2.303, WBIC_minus_FE0 = 2.038)
(k = 5, prob = 0.2007, FE_minus_FE0 = 2.396, BIC_minus_FE0 = 1.894, WBIC_minus_FE0 = 1.636)
(k = 6, prob = 0.1115, FE_minus_FE0 = 1.184, BIC_minus_FE0 = 0.681, WBIC_minus_FE0 = 0.416)
(k = 7, prob = 0.0425, FE_minus_FE0 = -0.865, BIC_minus_FE0 = -1.373, WBIC_minus_FE0 = -1.658)
(k = 8, prob = 0.0106, FE_minus_FE0 = -3.874, BIC_minus_FE0 = -4.393, WBIC_minus_FE0 = -4.723)
(k = 9, prob = 0.0016, FE_minus_FE0 = -8.154, BIC_minus_FE0 = -8.711, WBIC_minus_FE0 = -9.134)

Out[49]:

In [50]:

n = 20
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    F0 = FE0(q,n,k)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>x->rd(x,4), 
        FE_minus_FE0=FE(n,k)-F0|>rd, 
        BIC_minus_FE0=BIC(n,k)-F0|>rd,
        WBIC_minus_FE0=WBIC(n,k)-F0|>rd
    ) |> println
end

plot_BIC_WBIC_FE(n, w)

n = 20
w = 0.4
(k = 1, prob = 0.0005, FE_minus_FE0 = -9.764, BIC_minus_FE0 = -10.308, WBIC_minus_FE0 = -10.836)
(k = 2, prob = 0.0031, FE_minus_FE0 = -5.55, BIC_minus_FE0 = -6.056, WBIC_minus_FE0 = -6.464)
(k = 3, prob = 0.0123, FE_minus_FE0 = -2.469, BIC_minus_FE0 = -2.962, WBIC_minus_FE0 = -3.298)
(k = 4, prob = 0.035, FE_minus_FE0 = -0.179, BIC_minus_FE0 = -0.665, WBIC_minus_FE0 = -0.958)
(k = 5, prob = 0.0746, FE_minus_FE0 = 1.484, BIC_minus_FE0 = 1.001, WBIC_minus_FE0 = 0.736)
(k = 6, prob = 0.1244, FE_minus_FE0 = 2.611, BIC_minus_FE0 = 2.132, WBIC_minus_FE0 = 1.884)
(k = 7, prob = 0.1659, FE_minus_FE0 = 3.262, BIC_minus_FE0 = 2.784, WBIC_minus_FE0 = 2.549)
(k = 8, prob = 0.1797, FE_minus_FE0 = 3.473, BIC_minus_FE0 = 2.996, WBIC_minus_FE0 = 2.768)
(k = 9, prob = 0.1597, FE_minus_FE0 = 3.267, BIC_minus_FE0 = 2.79, WBIC_minus_FE0 = 2.567)
(k = 10, prob = 0.1171, FE_minus_FE0 = 2.656, BIC_minus_FE0 = 2.179, WBIC_minus_FE0 = 1.957)
(k = 11, prob = 0.071, FE_minus_FE0 = 1.645, BIC_minus_FE0 = 1.168, WBIC_minus_FE0 = 0.945)
(k = 12, prob = 0.0355, FE_minus_FE0 = 0.229, BIC_minus_FE0 = -0.248, WBIC_minus_FE0 = -0.476)
(k = 13, prob = 0.0146, FE_minus_FE0 = -1.603, BIC_minus_FE0 = -2.082, WBIC_minus_FE0 = -2.317)
(k = 14, prob = 0.0049, FE_minus_FE0 = -3.876, BIC_minus_FE0 = -4.356, WBIC_minus_FE0 = -4.603)
(k = 15, prob = 0.0013, FE_minus_FE0 = -6.626, BIC_minus_FE0 = -7.108, WBIC_minus_FE0 = -7.373)

Out[50]:

In [51]:

n = 50
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    F0 = FE0(q,n,k)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>x->rd(x,4), 
        FE_minus_FE0=FE(n,k)-F0|>rd, 
        BIC_minus_FE0=BIC(n,k)-F0|>rd,
        WBIC_minus_FE0=WBIC(n,k)-F0|>rd
    ) |> println
end

plot_BIC_WBIC_FE(n, w)

n = 50
w = 0.4
(k = 9, prob = 0.0005, FE_minus_FE0 = -6.863, BIC_minus_FE0 = -7.33, WBIC_minus_FE0 = -7.549)
(k = 10, prob = 0.0014, FE_minus_FE0 = -4.774, BIC_minus_FE0 = -5.24, WBIC_minus_FE0 = -5.447)
(k = 11, prob = 0.0035, FE_minus_FE0 = -2.935, BIC_minus_FE0 = -3.4, WBIC_minus_FE0 = -3.597)
(k = 12, prob = 0.0076, FE_minus_FE0 = -1.33, BIC_minus_FE0 = -1.794, WBIC_minus_FE0 = -1.983)
(k = 13, prob = 0.0147, FE_minus_FE0 = 0.057, BIC_minus_FE0 = -0.407, WBIC_minus_FE0 = -0.589)
(k = 14, prob = 0.026, FE_minus_FE0 = 1.235, BIC_minus_FE0 = 0.772, WBIC_minus_FE0 = 0.596)
(k = 15, prob = 0.0415, FE_minus_FE0 = 2.215, BIC_minus_FE0 = 1.752, WBIC_minus_FE0 = 1.581)
(k = 16, prob = 0.0606, FE_minus_FE0 = 3.004, BIC_minus_FE0 = 2.542, WBIC_minus_FE0 = 2.374)
(k = 17, prob = 0.0808, FE_minus_FE0 = 3.61, BIC_minus_FE0 = 3.147, WBIC_minus_FE0 = 2.983)
(k = 18, prob = 0.0987, FE_minus_FE0 = 4.037, BIC_minus_FE0 = 3.575, WBIC_minus_FE0 = 3.414)
(k = 19, prob = 0.1109, FE_minus_FE0 = 4.29, BIC_minus_FE0 = 3.828, WBIC_minus_FE0 = 3.67)
(k = 20, prob = 0.1146, FE_minus_FE0 = 4.374, BIC_minus_FE0 = 3.912, WBIC_minus_FE0 = 3.755)
(k = 21, prob = 0.1091, FE_minus_FE0 = 4.291, BIC_minus_FE0 = 3.829, WBIC_minus_FE0 = 3.674)
(k = 22, prob = 0.0959, FE_minus_FE0 = 4.044, BIC_minus_FE0 = 3.582, WBIC_minus_FE0 = 3.428)
(k = 23, prob = 0.0778, FE_minus_FE0 = 3.634, BIC_minus_FE0 = 3.172, WBIC_minus_FE0 = 3.019)
(k = 24, prob = 0.0584, FE_minus_FE0 = 3.063, BIC_minus_FE0 = 2.602, WBIC_minus_FE0 = 2.449)
(k = 25, prob = 0.0405, FE_minus_FE0 = 2.333, BIC_minus_FE0 = 1.871, WBIC_minus_FE0 = 1.718)
(k = 26, prob = 0.0259, FE_minus_FE0 = 1.442, BIC_minus_FE0 = 0.98, WBIC_minus_FE0 = 0.827)
(k = 27, prob = 0.0154, FE_minus_FE0 = 0.39, BIC_minus_FE0 = -0.071, WBIC_minus_FE0 = -0.224)
(k = 28, prob = 0.0084, FE_minus_FE0 = -0.822, BIC_minus_FE0 = -1.284, WBIC_minus_FE0 = -1.438)
(k = 29, prob = 0.0043, FE_minus_FE0 = -2.197, BIC_minus_FE0 = -2.658, WBIC_minus_FE0 = -2.814)
(k = 30, prob = 0.002, FE_minus_FE0 = -3.735, BIC_minus_FE0 = -4.197, WBIC_minus_FE0 = -4.354)
(k = 31, prob = 0.0009, FE_minus_FE0 = -5.441, BIC_minus_FE0 = -5.903, WBIC_minus_FE0 = -6.062)

Out[51]:

In [52]:

n = 100
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    F0 = FE0(q,n,k)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>x->rd(x,4), 
        FE_minus_FE0=FE(n,k)-F0|>rd, 
        BIC_minus_FE0=BIC(n,k)-F0|>rd,
        WBIC_minus_FE0=WBIC(n,k)-F0|>rd
    ) |> println
end

plot_BIC_WBIC_FE(n, w)

n = 100
w = 0.4
(k = 25, prob = 0.0006, FE_minus_FE0 = -4.908, BIC_minus_FE0 = -5.366, WBIC_minus_FE0 = -5.495)
(k = 26, prob = 0.0012, FE_minus_FE0 = -3.575, BIC_minus_FE0 = -4.033, WBIC_minus_FE0 = -4.16)
(k = 27, prob = 0.0022, FE_minus_FE0 = -2.346, BIC_minus_FE0 = -2.803, WBIC_minus_FE0 = -2.928)
(k = 28, prob = 0.0038, FE_minus_FE0 = -1.218, BIC_minus_FE0 = -1.675, WBIC_minus_FE0 = -1.798)
(k = 29, prob = 0.0063, FE_minus_FE0 = -0.189, BIC_minus_FE0 = -0.647, WBIC_minus_FE0 = -0.767)
(k = 30, prob = 0.01, FE_minus_FE0 = 0.742, BIC_minus_FE0 = 0.285, WBIC_minus_FE0 = 0.166)
(k = 31, prob = 0.0151, FE_minus_FE0 = 1.578, BIC_minus_FE0 = 1.121, WBIC_minus_FE0 = 1.004)
(k = 32, prob = 0.0217, FE_minus_FE0 = 2.321, BIC_minus_FE0 = 1.864, WBIC_minus_FE0 = 1.748)
(k = 33, prob = 0.0297, FE_minus_FE0 = 2.972, BIC_minus_FE0 = 2.515, WBIC_minus_FE0 = 2.4)
(k = 34, prob = 0.0391, FE_minus_FE0 = 3.532, BIC_minus_FE0 = 3.076, WBIC_minus_FE0 = 2.962)
(k = 35, prob = 0.0491, FE_minus_FE0 = 4.004, BIC_minus_FE0 = 3.547, WBIC_minus_FE0 = 3.434)
(k = 36, prob = 0.0591, FE_minus_FE0 = 4.387, BIC_minus_FE0 = 3.93, WBIC_minus_FE0 = 3.818)
(k = 37, prob = 0.0682, FE_minus_FE0 = 4.684, BIC_minus_FE0 = 4.227, WBIC_minus_FE0 = 4.116)
(k = 38, prob = 0.0754, FE_minus_FE0 = 4.894, BIC_minus_FE0 = 4.438, WBIC_minus_FE0 = 4.327)
(k = 39, prob = 0.0799, FE_minus_FE0 = 5.02, BIC_minus_FE0 = 4.563, WBIC_minus_FE0 = 4.454)
(k = 40, prob = 0.0812, FE_minus_FE0 = 5.062, BIC_minus_FE0 = 4.605, WBIC_minus_FE0 = 4.496)
(k = 41, prob = 0.0792, FE_minus_FE0 = 5.02, BIC_minus_FE0 = 4.564, WBIC_minus_FE0 = 4.455)
(k = 42, prob = 0.0742, FE_minus_FE0 = 4.896, BIC_minus_FE0 = 4.439, WBIC_minus_FE0 = 4.332)
(k = 43, prob = 0.0667, FE_minus_FE0 = 4.69, BIC_minus_FE0 = 4.233, WBIC_minus_FE0 = 4.126)
(k = 44, prob = 0.0576, FE_minus_FE0 = 4.402, BIC_minus_FE0 = 3.945, WBIC_minus_FE0 = 3.838)
(k = 45, prob = 0.0478, FE_minus_FE0 = 4.033, BIC_minus_FE0 = 3.576, WBIC_minus_FE0 = 3.469)
(k = 46, prob = 0.0381, FE_minus_FE0 = 3.583, BIC_minus_FE0 = 3.126, WBIC_minus_FE0 = 3.02)
(k = 47, prob = 0.0292, FE_minus_FE0 = 3.052, BIC_minus_FE0 = 2.596, WBIC_minus_FE0 = 2.489)
(k = 48, prob = 0.0215, FE_minus_FE0 = 2.441, BIC_minus_FE0 = 1.985, WBIC_minus_FE0 = 1.879)
(k = 49, prob = 0.0152, FE_minus_FE0 = 1.75, BIC_minus_FE0 = 1.294, WBIC_minus_FE0 = 1.188)
(k = 50, prob = 0.0103, FE_minus_FE0 = 0.98, BIC_minus_FE0 = 0.523, WBIC_minus_FE0 = 0.417)
(k = 51, prob = 0.0068, FE_minus_FE0 = 0.129, BIC_minus_FE0 = -0.328, WBIC_minus_FE0 = -0.434)
(k = 52, prob = 0.0042, FE_minus_FE0 = -0.802, BIC_minus_FE0 = -1.259, WBIC_minus_FE0 = -1.365)
(k = 53, prob = 0.0026, FE_minus_FE0 = -1.813, BIC_minus_FE0 = -2.27, WBIC_minus_FE0 = -2.376)
(k = 54, prob = 0.0015, FE_minus_FE0 = -2.905, BIC_minus_FE0 = -3.361, WBIC_minus_FE0 = -3.468)
(k = 55, prob = 0.0008, FE_minus_FE0 = -4.077, BIC_minus_FE0 = -4.533, WBIC_minus_FE0 = -4.64)

Out[52]:

In [53]:

n = 1000
w = 0.4
q = Bernoulli(w)
@show n
@show w
for k in sample_supp(q, n, 3)
    F0 = FE0(q,n,k)
    (
        k=k, 
        prob=pdf(Binomial(n, q.p), k)|>x->rd(x,4), 
        FE_minus_FE0=FE(n,k)-F0|>rd, 
        BIC_minus_FE0=BIC(n,k)-F0|>rd,
        WBIC_minus_FE0=WBIC(n,k)-F0|>rd
    ) |> println
end

plot_BIC_WBIC_FE(n, w; ylim=(-4, 8))

n = 1000
w = 0.4
(k = 353, prob = 0.0002, FE_minus_FE0 = -1.982, BIC_minus_FE0 = -2.434, WBIC_minus_FE0 = -2.459)
(k = 354, prob = 0.0003, FE_minus_FE0 = -1.585, BIC_minus_FE0 = -2.037, WBIC_minus_FE0 = -2.063)
(k = 355, prob = 0.0004, FE_minus_FE0 = -1.198, BIC_minus_FE0 = -1.65, WBIC_minus_FE0 = -1.675)
(k = 356, prob = 0.0004, FE_minus_FE0 = -0.819, BIC_minus_FE0 = -1.271, WBIC_minus_FE0 = -1.296)
(k = 357, prob = 0.0005, FE_minus_FE0 = -0.448, BIC_minus_FE0 = -0.9, WBIC_minus_FE0 = -0.926)
(k = 358, prob = 0.0006, FE_minus_FE0 = -0.087, BIC_minus_FE0 = -0.539, WBIC_minus_FE0 = -0.564)
(k = 359, prob = 0.0008, FE_minus_FE0 = 0.266, BIC_minus_FE0 = -0.186, WBIC_minus_FE0 = -0.211)
(k = 360, prob = 0.0009, FE_minus_FE0 = 0.61, BIC_minus_FE0 = 0.158, WBIC_minus_FE0 = 0.133)
(k = 361, prob = 0.0011, FE_minus_FE0 = 0.946, BIC_minus_FE0 = 0.493, WBIC_minus_FE0 = 0.468)
(k = 362, prob = 0.0013, FE_minus_FE0 = 1.272, BIC_minus_FE0 = 0.82, WBIC_minus_FE0 = 0.795)
(k = 363, prob = 0.0015, FE_minus_FE0 = 1.591, BIC_minus_FE0 = 1.138, WBIC_minus_FE0 = 1.113)
(k = 364, prob = 0.0017, FE_minus_FE0 = 1.9, BIC_minus_FE0 = 1.448, WBIC_minus_FE0 = 1.423)
(k = 365, prob = 0.002, FE_minus_FE0 = 2.201, BIC_minus_FE0 = 1.749, WBIC_minus_FE0 = 1.724)
(k = 366, prob = 0.0023, FE_minus_FE0 = 2.493, BIC_minus_FE0 = 2.041, WBIC_minus_FE0 = 2.016)
(k = 367, prob = 0.0026, FE_minus_FE0 = 2.777, BIC_minus_FE0 = 2.325, WBIC_minus_FE0 = 2.299)
(k = 368, prob = 0.003, FE_minus_FE0 = 3.052, BIC_minus_FE0 = 2.6, WBIC_minus_FE0 = 2.574)
(k = 369, prob = 0.0035, FE_minus_FE0 = 3.318, BIC_minus_FE0 = 2.866, WBIC_minus_FE0 = 2.841)
(k = 370, prob = 0.0039, FE_minus_FE0 = 3.576, BIC_minus_FE0 = 3.124, WBIC_minus_FE0 = 3.099)
(k = 371, prob = 0.0045, FE_minus_FE0 = 3.825, BIC_minus_FE0 = 3.373, WBIC_minus_FE0 = 3.348)
(k = 372, prob = 0.005, FE_minus_FE0 = 4.066, BIC_minus_FE0 = 3.614, WBIC_minus_FE0 = 3.589)
(k = 373, prob = 0.0056, FE_minus_FE0 = 4.298, BIC_minus_FE0 = 3.846, WBIC_minus_FE0 = 3.821)
(k = 374, prob = 0.0063, FE_minus_FE0 = 4.521, BIC_minus_FE0 = 4.069, WBIC_minus_FE0 = 4.044)
(k = 375, prob = 0.007, FE_minus_FE0 = 4.736, BIC_minus_FE0 = 4.284, WBIC_minus_FE0 = 4.259)
(k = 376, prob = 0.0078, FE_minus_FE0 = 4.943, BIC_minus_FE0 = 4.491, WBIC_minus_FE0 = 4.466)
(k = 377, prob = 0.0086, FE_minus_FE0 = 5.141, BIC_minus_FE0 = 4.689, WBIC_minus_FE0 = 4.664)
(k = 378, prob = 0.0094, FE_minus_FE0 = 5.33, BIC_minus_FE0 = 4.878, WBIC_minus_FE0 = 4.853)
(k = 379, prob = 0.0103, FE_minus_FE0 = 5.511, BIC_minus_FE0 = 5.059, WBIC_minus_FE0 = 5.034)
(k = 380, prob = 0.0112, FE_minus_FE0 = 5.683, BIC_minus_FE0 = 5.231, WBIC_minus_FE0 = 5.206)
(k = 381, prob = 0.0122, FE_minus_FE0 = 5.847, BIC_minus_FE0 = 5.395, WBIC_minus_FE0 = 5.37)
(k = 382, prob = 0.0132, FE_minus_FE0 = 6.003, BIC_minus_FE0 = 5.551, WBIC_minus_FE0 = 5.526)
(k = 383, prob = 0.0142, FE_minus_FE0 = 6.15, BIC_minus_FE0 = 5.698, WBIC_minus_FE0 = 5.673)
(k = 384, prob = 0.0152, FE_minus_FE0 = 6.288, BIC_minus_FE0 = 5.836, WBIC_minus_FE0 = 5.811)
(k = 385, prob = 0.0162, FE_minus_FE0 = 6.418, BIC_minus_FE0 = 5.966, WBIC_minus_FE0 = 5.941)
(k = 386, prob = 0.0172, FE_minus_FE0 = 6.54, BIC_minus_FE0 = 6.088, WBIC_minus_FE0 = 6.063)
(k = 387, prob = 0.0182, FE_minus_FE0 = 6.653, BIC_minus_FE0 = 6.201, WBIC_minus_FE0 = 6.176)
(k = 388, prob = 0.0192, FE_minus_FE0 = 6.758, BIC_minus_FE0 = 6.306, WBIC_minus_FE0 = 6.281)
(k = 389, prob = 0.0201, FE_minus_FE0 = 6.854, BIC_minus_FE0 = 6.402, WBIC_minus_FE0 = 6.377)
(k = 390, prob = 0.021, FE_minus_FE0 = 6.942, BIC_minus_FE0 = 6.49, WBIC_minus_FE0 = 6.465)
(k = 391, prob = 0.0218, FE_minus_FE0 = 7.021, BIC_minus_FE0 = 6.569, WBIC_minus_FE0 = 6.545)
(k = 392, prob = 0.0226, FE_minus_FE0 = 7.093, BIC_minus_FE0 = 6.64, WBIC_minus_FE0 = 6.616)
(k = 393, prob = 0.0233, FE_minus_FE0 = 7.155, BIC_minus_FE0 = 6.703, WBIC_minus_FE0 = 6.679)
(k = 394, prob = 0.0239, FE_minus_FE0 = 7.21, BIC_minus_FE0 = 6.758, WBIC_minus_FE0 = 6.733)
(k = 395, prob = 0.0245, FE_minus_FE0 = 7.256, BIC_minus_FE0 = 6.803, WBIC_minus_FE0 = 6.779)
(k = 396, prob = 0.0249, FE_minus_FE0 = 7.293, BIC_minus_FE0 = 6.841, WBIC_minus_FE0 = 6.816)
(k = 397, prob = 0.0253, FE_minus_FE0 = 7.322, BIC_minus_FE0 = 6.87, WBIC_minus_FE0 = 6.846)
(k = 398, prob = 0.0256, FE_minus_FE0 = 7.343, BIC_minus_FE0 = 6.891, WBIC_minus_FE0 = 6.867)
(k = 399, prob = 0.0257, FE_minus_FE0 = 7.356, BIC_minus_FE0 = 6.904, WBIC_minus_FE0 = 6.879)
(k = 400, prob = 0.0257, FE_minus_FE0 = 7.36, BIC_minus_FE0 = 6.908, WBIC_minus_FE0 = 6.883)
(k = 401, prob = 0.0257, FE_minus_FE0 = 7.356, BIC_minus_FE0 = 6.904, WBIC_minus_FE0 = 6.879)
(k = 402, prob = 0.0255, FE_minus_FE0 = 7.343, BIC_minus_FE0 = 6.891, WBIC_minus_FE0 = 6.867)
(k = 403, prob = 0.0252, FE_minus_FE0 = 7.322, BIC_minus_FE0 = 6.87, WBIC_minus_FE0 = 6.846)
(k = 404, prob = 0.0249, FE_minus_FE0 = 7.293, BIC_minus_FE0 = 6.841, WBIC_minus_FE0 = 6.817)
(k = 405, prob = 0.0244, FE_minus_FE0 = 7.256, BIC_minus_FE0 = 6.804, WBIC_minus_FE0 = 6.779)
(k = 406, prob = 0.0238, FE_minus_FE0 = 7.21, BIC_minus_FE0 = 6.758, WBIC_minus_FE0 = 6.734)
(k = 407, prob = 0.0232, FE_minus_FE0 = 7.156, BIC_minus_FE0 = 6.704, WBIC_minus_FE0 = 6.68)
(k = 408, prob = 0.0225, FE_minus_FE0 = 7.094, BIC_minus_FE0 = 6.642, WBIC_minus_FE0 = 6.617)
(k = 409, prob = 0.0217, FE_minus_FE0 = 7.023, BIC_minus_FE0 = 6.571, WBIC_minus_FE0 = 6.547)
(k = 410, prob = 0.0208, FE_minus_FE0 = 6.944, BIC_minus_FE0 = 6.492, WBIC_minus_FE0 = 6.468)
(k = 411, prob = 0.0199, FE_minus_FE0 = 6.857, BIC_minus_FE0 = 6.405, WBIC_minus_FE0 = 6.381)
(k = 412, prob = 0.019, FE_minus_FE0 = 6.762, BIC_minus_FE0 = 6.31, WBIC_minus_FE0 = 6.285)
(k = 413, prob = 0.018, FE_minus_FE0 = 6.658, BIC_minus_FE0 = 6.206, WBIC_minus_FE0 = 6.182)
(k = 414, prob = 0.017, FE_minus_FE0 = 6.546, BIC_minus_FE0 = 6.094, WBIC_minus_FE0 = 6.07)
(k = 415, prob = 0.016, FE_minus_FE0 = 6.426, BIC_minus_FE0 = 5.974, WBIC_minus_FE0 = 5.95)
(k = 416, prob = 0.015, FE_minus_FE0 = 6.298, BIC_minus_FE0 = 5.846, WBIC_minus_FE0 = 5.821)
(k = 417, prob = 0.014, FE_minus_FE0 = 6.161, BIC_minus_FE0 = 5.709, WBIC_minus_FE0 = 5.685)
(k = 418, prob = 0.0131, FE_minus_FE0 = 6.016, BIC_minus_FE0 = 5.564, WBIC_minus_FE0 = 5.54)
(k = 419, prob = 0.0121, FE_minus_FE0 = 5.863, BIC_minus_FE0 = 5.411, WBIC_minus_FE0 = 5.387)
(k = 420, prob = 0.0112, FE_minus_FE0 = 5.702, BIC_minus_FE0 = 5.25, WBIC_minus_FE0 = 5.226)
(k = 421, prob = 0.0102, FE_minus_FE0 = 5.532, BIC_minus_FE0 = 5.08, WBIC_minus_FE0 = 5.056)
(k = 422, prob = 0.0094, FE_minus_FE0 = 5.355, BIC_minus_FE0 = 4.903, WBIC_minus_FE0 = 4.878)
(k = 423, prob = 0.0085, FE_minus_FE0 = 5.169, BIC_minus_FE0 = 4.717, WBIC_minus_FE0 = 4.693)
(k = 424, prob = 0.0077, FE_minus_FE0 = 4.975, BIC_minus_FE0 = 4.523, WBIC_minus_FE0 = 4.498)
(k = 425, prob = 0.007, FE_minus_FE0 = 4.772, BIC_minus_FE0 = 4.32, WBIC_minus_FE0 = 4.296)
(k = 426, prob = 0.0063, FE_minus_FE0 = 4.562, BIC_minus_FE0 = 4.11, WBIC_minus_FE0 = 4.086)
(k = 427, prob = 0.0056, FE_minus_FE0 = 4.343, BIC_minus_FE0 = 3.891, WBIC_minus_FE0 = 3.867)
(k = 428, prob = 0.005, FE_minus_FE0 = 4.117, BIC_minus_FE0 = 3.664, WBIC_minus_FE0 = 3.64)
(k = 429, prob = 0.0045, FE_minus_FE0 = 3.882, BIC_minus_FE0 = 3.429, WBIC_minus_FE0 = 3.405)
(k = 430, prob = 0.004, FE_minus_FE0 = 3.638, BIC_minus_FE0 = 3.186, WBIC_minus_FE0 = 3.162)
(k = 431, prob = 0.0035, FE_minus_FE0 = 3.387, BIC_minus_FE0 = 2.935, WBIC_minus_FE0 = 2.911)
(k = 432, prob = 0.0031, FE_minus_FE0 = 3.128, BIC_minus_FE0 = 2.676, WBIC_minus_FE0 = 2.651)
(k = 433, prob = 0.0027, FE_minus_FE0 = 2.86, BIC_minus_FE0 = 2.408, WBIC_minus_FE0 = 2.384)
(k = 434, prob = 0.0023, FE_minus_FE0 = 2.584, BIC_minus_FE0 = 2.132, WBIC_minus_FE0 = 2.108)
(k = 435, prob = 0.002, FE_minus_FE0 = 2.3, BIC_minus_FE0 = 1.848, WBIC_minus_FE0 = 1.824)
(k = 436, prob = 0.0018, FE_minus_FE0 = 2.008, BIC_minus_FE0 = 1.556, WBIC_minus_FE0 = 1.532)
(k = 437, prob = 0.0015, FE_minus_FE0 = 1.708, BIC_minus_FE0 = 1.256, WBIC_minus_FE0 = 1.232)
(k = 438, prob = 0.0013, FE_minus_FE0 = 1.4, BIC_minus_FE0 = 0.948, WBIC_minus_FE0 = 0.924)
(k = 439, prob = 0.0011, FE_minus_FE0 = 1.083, BIC_minus_FE0 = 0.631, WBIC_minus_FE0 = 0.607)
(k = 440, prob = 0.0009, FE_minus_FE0 = 0.759, BIC_minus_FE0 = 0.307, WBIC_minus_FE0 = 0.283)
(k = 441, prob = 0.0008, FE_minus_FE0 = 0.426, BIC_minus_FE0 = -0.026, WBIC_minus_FE0 = -0.05)
(k = 442, prob = 0.0007, FE_minus_FE0 = 0.085, BIC_minus_FE0 = -0.367, WBIC_minus_FE0 = -0.391)
(k = 443, prob = 0.0006, FE_minus_FE0 = -0.263, BIC_minus_FE0 = -0.716, WBIC_minus_FE0 = -0.74)
(k = 444, prob = 0.0005, FE_minus_FE0 = -0.62, BIC_minus_FE0 = -1.073, WBIC_minus_FE0 = -1.097)
(k = 445, prob = 0.0004, FE_minus_FE0 = -0.986, BIC_minus_FE0 = -1.438, WBIC_minus_FE0 = -1.462)
(k = 446, prob = 0.0003, FE_minus_FE0 = -1.359, BIC_minus_FE0 = -1.811, WBIC_minus_FE0 = -1.835)
(k = 447, prob = 0.0003, FE_minus_FE0 = -1.74, BIC_minus_FE0 = -2.192, WBIC_minus_FE0 = -2.216)

Out[53]:

BICとWBICと自由エネルギーによって正しいモデル選択に失敗する確率¶

In [54]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_BIC = E(q, n, k -> BIC(n,k) < BIC0(q,n,k))
    (n=n, prob_fail_BIC=prob_fail_BIC|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_BIC = 0.101119)
(n = 32, prob_fail_BIC = 0.045704)
(n = 100, prob_fail_BIC = 0.031537)
(n = 316, prob_fail_BIC = 0.015748)
(n = 1000, prob_fail_BIC = 0.008909)
(n = 3162, prob_fail_BIC = 0.00437)
(n = 10000, prob_fail_BIC = 0.002434)
(n = 31623, prob_fail_BIC = 0.001308)
(n = 100000, prob_fail_BIC = 0.000685)
(n = 316228, prob_fail_BIC = 0.000375)
(n = 1000000, prob_fail_BIC = 0.000202)

In [55]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_WBIC = E(q, n, k -> WBIC(n,k) < WBIC0(q,n,k))
    (n=n, prob_fail_WBIC=prob_fail_WBIC|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_WBIC = 0.101119)
(n = 32, prob_fail_WBIC = 0.071091)
(n = 100, prob_fail_WBIC = 0.031537)
(n = 316, prob_fail_WBIC = 0.015748)
(n = 1000, prob_fail_WBIC = 0.008909)
(n = 3162, prob_fail_WBIC = 0.00437)
(n = 10000, prob_fail_WBIC = 0.002434)
(n = 31623, prob_fail_WBIC = 0.001308)
(n = 100000, prob_fail_WBIC = 0.000685)
(n = 316228, prob_fail_WBIC = 0.000375)
(n = 1000000, prob_fail_WBIC = 0.000202)

In [56]:

w = 0.4
@show w
q = Bernoulli(w)
for n in round.(Int, 10 .^ (1:0.5:6))
    prob_fail_FE = E(q, n, k -> FE(n,k) < FE0(q,n,k))
    (n=n, prob_fail_FE=prob_fail_FE|>x->rd(x,6)) |> println
end

w = 0.4
(n = 10, prob_fail_FE = 0.101119)
(n = 32, prob_fail_FE = 0.045704)
(n = 100, prob_fail_FE = 0.02478)
(n = 316, prob_fail_FE = 0.011404)
(n = 1000, prob_fail_FE = 0.00669)
(n = 3162, prob_fail_FE = 0.003469)
(n = 10000, prob_fail_FE = 0.001917)
(n = 31623, prob_fail_FE = 0.001007)
(n = 100000, prob_fail_FE = 0.00054)
(n = 316228, prob_fail_FE = 0.000292)
(n = 1000000, prob_fail_FE = 0.000159)

ベイズ統計と仮説検定や信頼区間との関係¶

In [57]:

# ジェネリックな信頼区間函数

is_in(x, CI) = CI[1] ≤ x ≤ CI[2]

function conf_int(pval_func, n, k; α=0.05)
    f(w) = pval_func(Bernoulli(w), n, k) - α
    CI = find_zeros(f, 0, 1)
    if k == 0
        [0.0, CI[1]]
    elseif k == n
        [CI[1], 1.0]
    else
        CI
    end
end

Out[57]:

conf_int (generic function with 1 method)

P値函数およびその類似物の定義¶

P値函数が定義されれば自動的に信頼区間も定義される.

In [58]:

# 二項分布の正規分布近似による仮説検定のP値

function pval_normal(q, n, k)
    q_sample = Binomial(n, q.p)
    μ = mean(q_sample)
    σ = std(q_sample)
    z = (k - μ)/σ
    2ccdf(Normal(), abs(z))
end

n = 10
[(k=k, conf_int=conf_int(pval_normal, n, k).|>rd) for k in 0:n]

Out[58]:

11-element Array{NamedTuple{(:k, :conf_int),Tuple{Int64,Array{Float64,1}}},1}:
 (k = 0, conf_int = [0.0, 0.278])  
 (k = 1, conf_int = [0.018, 0.404])
 (k = 2, conf_int = [0.057, 0.51]) 
 (k = 3, conf_int = [0.108, 0.603])
 (k = 4, conf_int = [0.168, 0.687])
 (k = 5, conf_int = [0.237, 0.763])
 (k = 6, conf_int = [0.313, 0.832])
 (k = 7, conf_int = [0.397, 0.892])
 (k = 8, conf_int = [0.49, 0.943]) 
 (k = 9, conf_int = [0.596, 0.982])
 (k = 10, conf_int = [0.722, 1.0])

In [59]:

# ベイズ信用区間

post_dist(n, k; a=0.5, b=0.5) = Beta(a+k, b+n-k)

cred_int(n, k; a=0.5, b=0.5, α=0.05) = quantile.(post_dist(n, k; a=a, b=b), [α/2, 1-α/2])

n = 10
[(k=k, cred_int=cred_int(n, k).|>rd) for k in 0:n]

Out[59]:

11-element Array{NamedTuple{(:k, :cred_int),Tuple{Int64,Array{Float64,1}}},1}:
 (k = 0, cred_int = [0.0, 0.217])  
 (k = 1, cred_int = [0.011, 0.381])
 (k = 2, cred_int = [0.044, 0.503])
 (k = 3, cred_int = [0.093, 0.606])
 (k = 4, cred_int = [0.153, 0.696])
 (k = 5, cred_int = [0.224, 0.776])
 (k = 6, cred_int = [0.304, 0.847])
 (k = 7, cred_int = [0.394, 0.907])
 (k = 8, cred_int = [0.497, 0.956])
 (k = 9, cred_int = [0.619, 0.989])
 (k = 10, cred_int = [0.783, 1.0])

In [60]:

# ベイズ信用区間の別の実装

# 以下の函数は事後分布を使ったP値の類似物
function pval_post(q, n, k; a=0.5, b=0.5)
    d = post_dist(n, k; a=a, b=b)
    min(2cdf(d, q.p), 2ccdf(d, q.p), 1)
end

# 上の計算結果と一致する
n = 10
[(k=k, conf_int=conf_int(pval_post, n, k).|>rd) for k in 0:n]

Out[60]:

11-element Array{NamedTuple{(:k, :conf_int),Tuple{Int64,Array{Float64,1}}},1}:
 (k = 0, conf_int = [0.0, 0.0])    
 (k = 1, conf_int = [0.011, 0.381])
 (k = 2, conf_int = [0.044, 0.503])
 (k = 3, conf_int = [0.093, 0.606])
 (k = 4, conf_int = [0.153, 0.696])
 (k = 5, conf_int = [0.224, 0.776])
 (k = 6, conf_int = [0.304, 0.847])
 (k = 7, conf_int = [0.394, 0.907])
 (k = 8, conf_int = [0.497, 0.956])
 (k = 9, conf_int = [0.619, 0.989])
 (k = 10, conf_int = [0.783, 1.0])

In [61]:

# AICを使ったP値の類似物 (対数尤度比のχ²検定のP値)

function pval_AIC(q, n, k)
    x = AIC0(q, n, k) - (AIC(n, k) - 2)
    ccdf(Chisq(1), x)
end

n = 10
[(k=k, conf_int=conf_int(pval_AIC, n, k).|>rd) for k in 0:n]

Out[61]:

11-element Array{NamedTuple{(:k, :conf_int),Tuple{Int64,Array{Float64,1}}},1}:
 (k = 0, conf_int = [0.0, 0.175])  
 (k = 1, conf_int = [0.006, 0.372])
 (k = 2, conf_int = [0.036, 0.499])
 (k = 3, conf_int = [0.085, 0.607])
 (k = 4, conf_int = [0.146, 0.7])  
 (k = 5, conf_int = [0.218, 0.782])
 (k = 6, conf_int = [0.3, 0.854])  
 (k = 7, conf_int = [0.393, 0.915])
 (k = 8, conf_int = [0.501, 0.964])
 (k = 9, conf_int = [0.628, 0.994])
 (k = 10, conf_int = [0.825, 1.0])

In [62]:

# WAICを使ったP値の類似物

function pval_WAIC(q, n, k; a=0.5, b=0.5)
    #x = WAIC0(q, n, k; a=a, b=b) - (WAIC(n, k; a=a, b=b) - 2)
    x = WAIC0(q, n, k; a=a, b=b) - 2n*T(n, k; a=a, b=b)
    ccdf(Chisq(1), x)
end

n = 10
[(k=k, conf_int=conf_int(pval_WAIC, n, k).|>rd) for k in 0:n]

Out[62]:

11-element Array{NamedTuple{(:k, :conf_int),Tuple{Int64,Array{Float64,1}}},1}:
 (k = 0, conf_int = [0.0, 0.212])  
 (k = 1, conf_int = [0.006, 0.377])
 (k = 2, conf_int = [0.036, 0.501])
 (k = 3, conf_int = [0.084, 0.607])
 (k = 4, conf_int = [0.146, 0.7])  
 (k = 5, conf_int = [0.218, 0.782])
 (k = 6, conf_int = [0.3, 0.854])  
 (k = 7, conf_int = [0.393, 0.916])
 (k = 8, conf_int = [0.499, 0.964])
 (k = 9, conf_int = [0.623, 0.994])
 (k = 10, conf_int = [0.788, 1.0])

In [63]:

n = 40
w = 0.29
q = Bernoulli(w)
sample = rand(Binomial(n, q.p), 10^5)
y1 = WAIC0.(q, n, sample) - WAIC.(n, sample) .+ 2
y2 = WAIC0.(q, n, sample) - 2n*T.(n, sample)
P1 = histogram(y1, xlim=(-1, 10), bin=-1:1:10, alpha=0.5, norm=true, label="WAIC0 - (WAIC - 2)")
plot!(x->pdf(Chisq(1),x), 0.2, 10; label="Chisq(1)")
P2 = histogram(y2, xlim=(-1, 10), bin=-1:1:10, alpha=0.5, norm=true, label="WAIC0 - 2n T")
plot!(x->pdf(Chisq(1),x), 0.2, 10; label="Chisq(1)")
plot(P1, P2, size=(800, 250))

Out[63]:

In [64]:

n = 100
w = 0.3
q = Bernoulli(w)
@show n
@show w
[(k=k, pval_normal=pval_normal(q, n, k)|>rd, pval_post=pval_post(q, n, k)|>rd, pval_AIC=pval_AIC(q, n, k)|>rd, pval_WAIC=pval_WAIC(q, n, k)|>rd) for k in sample_supp(q, n, 3)]

n = 100
w = 0.3

Out[64]:

29-element Array{NamedTuple{(:k, :pval_normal, :pval_post, :pval_AIC, :pval_WAIC),Tuple{Int64,Float64,Float64,Float64,Float64}},1}:
 (k = 16, pval_normal = 0.002, pval_post = 0.001, pval_AIC = 0.001, pval_WAIC = 0.001)
 (k = 17, pval_normal = 0.005, pval_post = 0.003, pval_AIC = 0.003, pval_WAIC = 0.003)
 (k = 18, pval_normal = 0.009, pval_post = 0.006, pval_AIC = 0.006, pval_WAIC = 0.006)
 (k = 19, pval_normal = 0.016, pval_post = 0.013, pval_AIC = 0.012, pval_WAIC = 0.012)
 (k = 20, pval_normal = 0.029, pval_post = 0.024, pval_AIC = 0.023, pval_WAIC = 0.023)
 (k = 21, pval_normal = 0.05, pval_post = 0.044, pval_AIC = 0.042, pval_WAIC = 0.042) 
 (k = 22, pval_normal = 0.081, pval_post = 0.075, pval_AIC = 0.072, pval_WAIC = 0.072)
 (k = 23, pval_normal = 0.127, pval_post = 0.121, pval_AIC = 0.117, pval_WAIC = 0.117)
 (k = 24, pval_normal = 0.19, pval_post = 0.186, pval_AIC = 0.181, pval_WAIC = 0.181) 
 (k = 25, pval_normal = 0.275, pval_post = 0.274, pval_AIC = 0.267, pval_WAIC = 0.267)
 (k = 26, pval_normal = 0.383, pval_post = 0.385, pval_AIC = 0.376, pval_WAIC = 0.377)
 (k = 27, pval_normal = 0.513, pval_post = 0.518, pval_AIC = 0.508, pval_WAIC = 0.51) 
 (k = 28, pval_normal = 0.663, pval_post = 0.671, pval_AIC = 0.66, pval_WAIC = 0.662) 
 (k = 29, pval_normal = 0.827, pval_post = 0.838, pval_AIC = 0.827, pval_WAIC = 0.83) 
 (k = 30, pval_normal = 1.0, pval_post = 0.988, pval_AIC = 1.0, pval_WAIC = 1.0)      
 (k = 31, pval_normal = 0.827, pval_post = 0.817, pval_AIC = 0.828, pval_WAIC = 0.831)
 (k = 32, pval_normal = 0.663, pval_post = 0.654, pval_AIC = 0.664, pval_WAIC = 0.666)
 (k = 33, pval_normal = 0.513, pval_post = 0.508, pval_AIC = 0.516, pval_WAIC = 0.517)
 (k = 34, pval_normal = 0.383, pval_post = 0.381, pval_AIC = 0.388, pval_WAIC = 0.389)
 (k = 35, pval_normal = 0.275, pval_post = 0.276, pval_AIC = 0.282, pval_WAIC = 0.282)
 (k = 36, pval_normal = 0.19, pval_post = 0.193, pval_AIC = 0.198, pval_WAIC = 0.198) 
 (k = 37, pval_normal = 0.127, pval_post = 0.131, pval_AIC = 0.134, pval_WAIC = 0.134)
 (k = 38, pval_normal = 0.081, pval_post = 0.085, pval_AIC = 0.088, pval_WAIC = 0.088)
 (k = 39, pval_normal = 0.05, pval_post = 0.054, pval_AIC = 0.055, pval_WAIC = 0.055) 
 (k = 40, pval_normal = 0.029, pval_post = 0.033, pval_AIC = 0.034, pval_WAIC = 0.034)
 (k = 41, pval_normal = 0.016, pval_post = 0.019, pval_AIC = 0.02, pval_WAIC = 0.02)  
 (k = 42, pval_normal = 0.009, pval_post = 0.011, pval_AIC = 0.011, pval_WAIC = 0.011)
 (k = 43, pval_normal = 0.005, pval_post = 0.006, pval_AIC = 0.006, pval_WAIC = 0.006)
 (k = 44, pval_normal = 0.002, pval_post = 0.003, pval_AIC = 0.003, pval_WAIC = 0.003)

In [65]:

n = 30
[
    (
        k=k, 
        Normal=conf_int(pval_normal, n, k).|>rd, 
        Bayes=cred_int(n, k).|>rd, 
        AIC=conf_int(pval_AIC, n, k).|>rd, 
        WAIC=conf_int(pval_WAIC, n, k).|>rd
    ) for k in 0:n
]

Out[65]:

31-element Array{NamedTuple{(:k, :Normal, :Bayes, :AIC, :WAIC),Tuple{Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1}}},1}:
 (k = 0, Normal = [0.0, 0.114], Bayes = [0.0, 0.08], AIC = [0.0, 0.062], WAIC = [0.0, 0.077])          
 (k = 1, Normal = [0.006, 0.167], Bayes = [0.004, 0.145], AIC = [0.002, 0.139], WAIC = [0.002, 0.142]) 
 (k = 2, Normal = [0.018, 0.213], Bayes = [0.014, 0.197], AIC = [0.011, 0.192], WAIC = [0.011, 0.194]) 
 (k = 3, Normal = [0.035, 0.256], Bayes = [0.029, 0.243], AIC = [0.026, 0.239], WAIC = [0.026, 0.24])  
 (k = 4, Normal = [0.053, 0.297], Bayes = [0.047, 0.287], AIC = [0.043, 0.283], WAIC = [0.043, 0.284]) 
 (k = 5, Normal = [0.073, 0.336], Bayes = [0.067, 0.327], AIC = [0.063, 0.325], WAIC = [0.063, 0.325]) 
 (k = 6, Normal = [0.095, 0.373], Bayes = [0.088, 0.367], AIC = [0.085, 0.364], WAIC = [0.085, 0.365]) 
 (k = 7, Normal = [0.118, 0.409], Bayes = [0.111, 0.404], AIC = [0.108, 0.403], WAIC = [0.108, 0.403]) 
 (k = 8, Normal = [0.142, 0.444], Bayes = [0.135, 0.441], AIC = [0.132, 0.44], WAIC = [0.132, 0.44])   
 (k = 9, Normal = [0.167, 0.479], Bayes = [0.16, 0.477], AIC = [0.157, 0.476], WAIC = [0.157, 0.476])  
 (k = 10, Normal = [0.192, 0.512], Bayes = [0.186, 0.511], AIC = [0.183, 0.511], WAIC = [0.183, 0.511])
 (k = 11, Normal = [0.219, 0.545], Bayes = [0.213, 0.545], AIC = [0.21, 0.545], WAIC = [0.21, 0.545])  
 (k = 12, Normal = [0.246, 0.577], Bayes = [0.24, 0.578], AIC = [0.238, 0.578], WAIC = [0.238, 0.578]) 
 (k = 13, Normal = [0.274, 0.608], Bayes = [0.269, 0.61], AIC = [0.267, 0.611], WAIC = [0.267, 0.611]) 
 (k = 14, Normal = [0.302, 0.639], Bayes = [0.298, 0.641], AIC = [0.296, 0.642], WAIC = [0.296, 0.642])
 (k = 15, Normal = [0.332, 0.668], Bayes = [0.328, 0.672], AIC = [0.327, 0.673], WAIC = [0.327, 0.673])
 (k = 16, Normal = [0.361, 0.698], Bayes = [0.359, 0.702], AIC = [0.358, 0.704], WAIC = [0.358, 0.704])
 (k = 17, Normal = [0.392, 0.726], Bayes = [0.39, 0.731], AIC = [0.389, 0.733], WAIC = [0.389, 0.733]) 
 (k = 18, Normal = [0.423, 0.754], Bayes = [0.422, 0.76], AIC = [0.422, 0.762], WAIC = [0.422, 0.762]) 
 (k = 19, Normal = [0.455, 0.781], Bayes = [0.455, 0.787], AIC = [0.455, 0.79], WAIC = [0.455, 0.79])  
 (k = 20, Normal = [0.488, 0.808], Bayes = [0.489, 0.814], AIC = [0.489, 0.817], WAIC = [0.489, 0.817])
 (k = 21, Normal = [0.521, 0.833], Bayes = [0.523, 0.84], AIC = [0.524, 0.843], WAIC = [0.524, 0.843]) 
 (k = 22, Normal = [0.556, 0.858], Bayes = [0.559, 0.865], AIC = [0.56, 0.868], WAIC = [0.56, 0.868])  
 (k = 23, Normal = [0.591, 0.882], Bayes = [0.596, 0.889], AIC = [0.597, 0.892], WAIC = [0.597, 0.892])
 (k = 24, Normal = [0.627, 0.905], Bayes = [0.633, 0.912], AIC = [0.636, 0.915], WAIC = [0.635, 0.915])
 (k = 25, Normal = [0.664, 0.927], Bayes = [0.673, 0.933], AIC = [0.675, 0.937], WAIC = [0.675, 0.937])
 (k = 26, Normal = [0.703, 0.947], Bayes = [0.713, 0.953], AIC = [0.717, 0.957], WAIC = [0.716, 0.957])
 (k = 27, Normal = [0.744, 0.965], Bayes = [0.757, 0.971], AIC = [0.761, 0.974], WAIC = [0.76, 0.974]) 
 (k = 28, Normal = [0.787, 0.982], Bayes = [0.803, 0.986], AIC = [0.808, 0.989], WAIC = [0.806, 0.989])
 (k = 29, Normal = [0.833, 0.994], Bayes = [0.855, 0.996], AIC = [0.861, 0.998], WAIC = [0.858, 0.998])
 (k = 30, Normal = [0.886, 1.0], Bayes = [0.92, 1.0], AIC = [0.938, 1.0], WAIC = [0.923, 1.0])

P値函数のアニメーション¶

In [66]:

function animate_pvalue_fuctions(n; fps=n/10, f=fill(true,4), pal=palette(:default))
    anim = @animate for k in 0:n
        w = 0.01:0.01:0.99
        y_normal = pval_normal.(Bernoulli.(w), n, k)
        y_Bayes  = pval_post.(Bernoulli.(w), n, k)
        y_AIC    = pval_AIC.(Bernoulli.(w), n, k)
        y_WAIC   = pval_WAIC.(Bernoulli.(w), n, k)
        P = plot(title="p-value functions for n=$n, k=" * @sprintf("%03d",k), titlefontsize=10)
        plot!(ylim=(-0.02, 1.02))
        plot!(xtick=0:0.1:1, ytick=0:0.1:1)
        f[1] && plot!(w, y_normal; label="Normal", ls=:solid, color=pal[1])
        f[2] && plot!(w, y_Bayes; label="Posterior", ls=:dash, color=pal[2])
        f[3] && plot!(w, y_AIC; label="AIC", ls=:dashdot, color=pal[3])
        f[4] && plot!(w, y_WAIC; label="WAIC", ls=:dot, lw=1.3, color=pal[4])
        plot(P; size=(500, 350))
    end
    pyplotclf()

    g = Int.(f)
    gifname = "pvalue_functions_$(n)_$(g[1])$(g[2])$(g[3])$(g[4]).gif"
    gif(anim, gifname; fps=fps)
    displayfile("image/gif", gifname)
end

Out[66]:

animate_pvalue_fuctions (generic function with 1 method)

In [67]:

@time animate_pvalue_fuctions(10)

  7.061318 seconds (11.13 M allocations: 664.561 MiB, 3.31% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_10_1111.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [68]:

@time animate_pvalue_fuctions(20)

  4.337529 seconds (423.63 k allocations: 24.159 MiB, 0.48% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_20_1111.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [69]:

@time animate_pvalue_fuctions(50)

 10.304655 seconds (1.02 M allocations: 57.770 MiB, 0.26% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_50_1111.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [70]:

@time animate_pvalue_fuctions(100)

 20.655136 seconds (2.03 M allocations: 114.074 MiB, 0.24% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_100_1111.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [71]:

@time animate_pvalue_fuctions(10; f=[true, true, false, false])

  2.250397 seconds (638.12 k allocations: 32.414 MiB)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_10_1100.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [72]:

@time animate_pvalue_fuctions(20; f=[true, true, false, false])

  3.912727 seconds (345.49 k allocations: 19.114 MiB, 0.53% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_20_1100.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [73]:

@time animate_pvalue_fuctions(50; f=[true, true, false, false])

  9.416736 seconds (831.77 k allocations: 45.685 MiB)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_50_1100.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [74]:

@time animate_pvalue_fuctions(100; f=[true, true, false, false])

 18.367178 seconds (1.64 M allocations: 90.194 MiB, 0.27% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0034\pvalue_functions_100_1100.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\WwFyB\src\animation.jl:98

In [ ]: