黒木玄
2017-11-16
このノートブックは
http://nbviewer.jupyter.org/gist/genkuroki/d688c6e5e8c427dba9856d37a12bbee1
で公開されているJupyter notebookで作成したデータファイルを解析するためのノートブックである.
使い方
(1) https://genkuroki.github.io/documents/Jupyter/#2017-11-16 から次の3つのファイルをダウンロードしてこのノートブックと同じ場所に置いておく.
(2) このノートブックを実行する.
n = 32
32
@time begin
using Mamba
using KernelDensity
function makefunc_pdfkde(X)
ik = InterpKDE(kde(X))
pdfkde(x) = pdf(ik, x)
return pdfkde
end
function makefunc_pdfkde(X,Y)
ik = InterpKDE(kde((X,Y)))
pdfkde(x, y) = pdf(ik, x, y)
return pdfkde
end
using Optim
optim_options = Optim.Options(
store_trace = true,
extended_trace = true
)
using QuadGK
import PyPlot
plt = PyPlot
using Distributions
@everywhere GTDist(μ, ρ, ν) = LocationScale(Float64(μ), Float64(ρ), TDist(Float64(ν)))
@everywhere GTDist(ρ, ν) = GTDist(zero(ρ), ρ, ν)
using JLD2
using FileIO
end
macro sum(f_k, k, itr)
quote
begin
s = zero(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))
for $(esc(k)) in $(esc(itr))
s += $(esc(f_k))
end
s
end
end
end
18.167492 seconds (14.84 M allocations: 958.267 MiB, 3.11% gc time)
@sum (macro with 1 method)
s(symbol, suffix) = eval(Symbol(symbol, :_, suffix))
function plotTh12(modelname::String; bins=25)
GL = s(:GL, modelname)
T = s(:T, modelname)
V = s(:V, modelname)
WAIC = s(:WAIC, modelname)
LOOCV = s(:LOOCV, modelname)
WBIC = s(:WBIC, modelname)
FreeEnergy = s(:FreeEnergy, modelname)
WG = WAIC - GL
LG = LOOCV - GL
LW = LOOCV - WAIC
KL = GL - Shannon
TT = T - T_true
WT = WAIC - T_true
KW = KL + WT
println("\n========== $modelname model (n = $n)")
println("2λ (std) = $(mean(KW)/2) ($(std(KW)/2))")
println("2ν (std) = $(mean(V)/2) ($(std(V)/2))")
println("="^70)
println("Name | mean (std)")
println("-"^70)
println("T | $(mean(T)) ($(std(T)))")
println("V | $(mean(V)) ($(std(V)))")
println("WAIC = T + V | $(mean(WAIC)) ($(std(WAIC)))")
println("LOOCV | $(mean(LOOCV)) ($(std(LOOCV)))")
println("GL | $(mean(GL)) ($(std(GL)))")
println("-"^70)
println("WAIC - GL | $(mean(WG)) ($(std(WG)))")
println("LOOCV - GL | $(mean(LG)) ($(std(LG)))")
println("LOOCV - WAIC | $(mean(LW)) ($(std(LW)))")
println("-"^70)
println("KL = GL - S | $(mean(KL)) ($(std(KL)))")
println("TT = T - T_true | $(mean(TT)) ($(std(TT)))")
println("WT = WAIC - T_true | $(mean(WT)) ($(std(WT)))")
println("KW = KL + WT | $(mean(KW)) ($(std(KW)))")
println("-"^70)
println("WBIC | $(mean(WBIC)) ($(std(WBIC)))")
println("FreeEnergy | $(mean(FreeEnergy)) ($(std(FreeEnergy)))")
println("WBIC - FreeEnergy | $(mean(WBIC-FreeEnergy)) ($(std(WBIC-FreeEnergy)))")
println("WBIC - T_true | $(mean(WBIC-T_true)) ($(std(WBIC-T_true)))")
println("FreeEnergy - T_true | $(mean(FreeEnergy-T_true)) ($(std(FreeEnergy-T_true)))")
println("="^70)
println("correlation of TT and KL = $(cor(TT, KL))")
println("correlation of WT and KL = $(cor(WT, KL))")
println("="^70)
println()
sleep(0.1)
plt.figure(figsize=(8,8))
plt.subplot(321)
x = linspace(minimum(WAIC), maximum(WAIC), 200)
plt.plt[:hist](WAIC, normed=true, bins=bins, alpha=0.5, label="WAIC")
plt.plot(x, pdf.(Chisq(n), x-(Shannon-n+mean(KW))), label="χ²($n) + const")
plt.grid(ls=":")
plt.legend(fontsize=8)
plt.title("WAIC of $modelname, n=$n, Nsims=$Nsims", fontsize=10)
plt.subplot(322)
x = linspace(minimum(T_true), maximum(T_true), 200)
plt.plt[:hist](T_true, normed=true, bins=bins, alpha=0.5, label="T_true")
plt.plot(x, pdf.(Chisq(n), x - (Shannon-n)), label="χ²($n) + const")
plt.grid(ls=":")
plt.legend(fontsize=8)
plt.title("T_true of the true dist, n=$n, Nsims=$Nsims", fontsize=10)
plt.subplot(323)
plt.plt[:hist](KL, normed=true, bins=bins, alpha=0.5, label="KL = GL \$-\$ S")
plt.grid(ls=":")
plt.legend()
plt.title("KL = GL \$-\$ Shannon for $modelname, n=$n", fontsize=10)
plt.subplot(324)
plt.plt[:hist](TT, normed=true, bins=bins, alpha=0.5, label="TT = T \$-\$ T_true")
plt.grid(ls=":")
plt.legend()
plt.title("TT = T \$-\$ T_true for $modelname, n=$n", fontsize=10)
plt.subplot(325)
plt.plt[:hist](WT, normed=true, bins=bins, alpha=0.5, label="WT = WAIC \$-\$ T_true")
plt.grid(ls=":")
plt.legend()
plt.title("WT = WAIC \$-\$ T_true for $modelname, n=$n", fontsize=10)
plt.subplot(326)
plt.plt[:hist](KW, normed=true, bins=bins, alpha=0.5, label="KW = KL + WT")
plt.axvline(mean(KW), color="k", ls=":", label="mean")
plt.grid(ls=":")
plt.legend(fontsize=8)
plt.title("KW = KL + WT for $modelname, n=$n", fontsize=10)
plt.tight_layout()
end
function plotModelSelection(criterion::String, leftmodel::String, rightmodel::String;
title="Model Selection by $criterion (n=$n)",
xlabel="$leftmodel ←―――――――――――――――――→ $rightmodel")
IC_left = s(criterion, leftmodel)
IC_right = s(criterion, rightmodel)
Diff = IC_left - IC_right
Percent_left = 100*count(x -> x < 0.0, Diff)/length(Diff)
Percent_right = 100*count(x -> x ≥ 0.0, Diff)/length(Diff)
rmin = minDiff = minimum(Diff)
rmax = maxDiff = maximum(Diff)
(rmin > -0.3*rmax) && (rmin = -0.3*rmax)
(rmax < -0.3*rmin) && (rmax = -0.3*rmin)
h = plt.plt[:hist](Diff, normed=true, bins=50, alpha=0.5, range=(rmin,rmax), label="diff of $criterion")
plt.xlabel(xlabel)
plt.ylabel("probability density")
plt.axvline(0, color="red", ls="--")
xann1 = rmin
yann1 = 0.9*maximum(h[1])
plt.annotate("$Percent_left%",
xytext=(0.9*xann1, 0.6*yann1),
xy=(0.4*minDiff, 0.2*yann1),
arrowprops=Dict("arrowstyle" => "->", "color" => "red"),
color="red")
xann2 = 0.7*rmax
yann2 = 0.9*maximum(h[1])
plt.annotate("$Percent_right%",
xytext=(0.9*xann2, 0.6*yann2),
xy=(0.4*maxDiff, 0.2*yann2),
arrowprops=Dict("arrowstyle" => "->", "color" => "red"),
color="red")
plt.legend()
plt.grid(ls=":")
plt.title(title, fontsize=10)
end
function plotAllModelSelections(criterion::String)
plt.figure(figsize=(10,7))
plt.subplot(221); plotModelSelection(criterion, "normal1", "mixnormal")
plt.subplot(222); plotModelSelection(criterion, "normal", "mixnormal")
plt.subplot(223); plotModelSelection(criterion, "normal", "normal1")
plt.tight_layout()
end
plotAllModelSelections (generic function with 1 method)
#data = load("DataStandardNormalSample8Case4649.jld2")
#data = load("DataStandardNormalSample32Case4649.jld2")
#data = load("DataStandardNormalSample128Case4649.jld2")
data = load("DataStandardNormalSample$(n)Case4649.jld2")
for v in sort(collect(keys(data)))
ex = parse("$v = data[\"$v\"]")
#println(ex)
eval(ex)
end
plotTh12("mixnormal")
plotTh12("normal1")
plotTh12("normal")
plotAllModelSelections("T")
plotAllModelSelections("GL")
plotAllModelSelections("WAIC")
plotAllModelSelections("LOOCV")
plotAllModelSelections("WBIC")
plotAllModelSelections("FreeEnergy");
========== mixnormal model (n = 32) 2λ (std) = 1.155175814734873 (0.1813903105053722) 2ν (std) = 1.0742079398050077 (0.4299022217524213) ====================================================================== Name | mean (std) ---------------------------------------------------------------------- T | 89.71453458532856 (6.604903567637937) V | 2.1484158796100155 (0.8598044435048426) WAIC = T + V | 91.86295046493856 (7.371807010659133) LOOCV | 91.86165485438019 (7.365566321980648) GL | 91.97317169532286 (1.2051202161380583) ---------------------------------------------------------------------- WAIC - GL | -0.11022123038427012 (7.245313316555237) LOOCV - GL | -0.11151684094266495 (7.239324044442848) LOOCV - WAIC | -0.0012956105583948555 (0.02973944581736434) ---------------------------------------------------------------------- KL = GL - S | 1.161105570223802 (1.2051202161380585) TT = T - T_true | -0.9991698203640712 (1.8341441540735908) WT = WAIC - T_true | 1.1492460592459444 (1.5218424853234411) KW = KL + WT | 2.310351629469746 (0.3627806210107444) ---------------------------------------------------------------------- WBIC | 93.12501637677528 (6.984221756458585) FreeEnergy | 92.40576360391438 (6.962418987983377) WBIC - FreeEnergy | 0.7192527728609136 (0.4183977423682536) WBIC - T_true | 2.411311971082669 (1.991077828522291) FreeEnergy - T_true | 1.692059198221755 (1.7877247594393728) ====================================================================== correlation of TT and KL = -0.9136491693919635 correlation of WT and KL = -0.9914676309931532 ====================================================================== ========== normal1 model (n = 32) 2λ (std) = 0.9330719597667703 (0.15076764029777262) 2ν (std) = 0.9482629541999565 (0.23651080855026488) ====================================================================== Name | mean (std) ---------------------------------------------------------------------- T | 89.71321551952528 (7.368064527839469) V | 1.896525908399913 (0.47302161710052976) WAIC = T + V | 91.60974142792519 (7.829031589695603) LOOCV | 91.61142810745376 (7.829452021917573) GL | 91.78217302240004 (1.3399528600180424) ---------------------------------------------------------------------- WAIC - GL | -0.17243159447482673 (7.9773306007179405) LOOCV - GL | -0.17074491494625815 (7.977763591516943) LOOCV - WAIC | 0.0016866795285685186 (0.0033420913200539167) ---------------------------------------------------------------------- KL = GL - S | 0.9701068973009779 (1.3399528600180424) TT = T - T_true | -1.0004888861673502 (1.4823572048977707) WT = WAIC - T_true | 0.8960370222325631 (1.5004116791513726) KW = KL + WT | 1.8661439195335405 (0.30153528059554524) ---------------------------------------------------------------------- WBIC | 92.57066964653924 (7.766540471455818) FreeEnergy | 91.97466013585884 (7.666203176930305) WBIC - FreeEnergy | 0.5960095106804452 (0.5864672756101755) WBIC - T_true | 1.8569652408466168 (2.0681075319960573) FreeEnergy - T_true | 1.2609557301661718 (1.7162923302533113) ====================================================================== correlation of TT and KL = -0.9876854839085023 correlation of WT and KL = -0.9837908203615242 ====================================================================== ========== normal model (n = 32) 2λ (std) = 1.9557216808208768 (0.41838333719197185) 2ν (std) = 1.80091166751125 (0.33102869477470764) ====================================================================== Name | mean (std) ---------------------------------------------------------------------- T | 89.06190569821592 (8.108806933610438) V | 3.6018233350225 (0.6620573895494153) WAIC = T + V | 92.66372903323841 (8.116944745587524) LOOCV | 92.71657877982683 (8.122926059015976) GL | 92.77348485919502 (1.9653333091267573) ---------------------------------------------------------------------- WAIC - GL | -0.10975582595658047 (8.558186765597165) LOOCV - GL | -0.05690607936820591 (8.564761570196836) LOOCV - WAIC | 0.05284974658837476 (0.08431686367806858) ---------------------------------------------------------------------- KL = GL - S | 1.961418734095962 (1.9653333091267577) TT = T - T_true | -1.6517987074767075 (2.083601038695401) WT = WAIC - T_true | 1.950024627545792 (2.1290880841650597) KW = KL + WT | 3.9114433616417537 (0.8367666743839437) ---------------------------------------------------------------------- WBIC | 94.91599743621809 (8.245736767533247) FreeEnergy | 93.50820328199669 (8.10096824594339) WBIC - FreeEnergy | 1.4077941542214105 (1.0038081392697724) WBIC - T_true | 4.202293030525466 (3.052600632531801) FreeEnergy - T_true | 2.7944988763040546 (2.4243462877296116) ====================================================================== correlation of TT and KL = -0.9351331311893791 correlation of WT and KL = -0.9195382809923446 ======================================================================
function plotAllTests(criterion)
model = ["normal1", "mixnormal", "normal"]
plt.figure(figsize=(10,7))
m = 0
for i in 1:3
for j in i+1:3
m += 1
Left = eval(Symbol(criterion, :_, model[i]))
Right = eval(Symbol(criterion, :_, model[j]))
Diff = Left - Right
a_Right = quantile(Diff, 0.05)
a_Left = quantile(Diff, 0.95)
sleep(0.1)
plt.subplot(2,2,m)
plt.plt[:hist](Diff, normed=true, bins=30, alpha=0.5)
plt.axvline(0.0, color="red", ls="--")
plt.axvline(a_Right, color="red", ls=":")
plt.axvline(a_Left, color="red", ls=":")
plt.xlabel("$(model[i]) ←――――――――――――――――→ $(model[j])")
plt.ylabel("probability density")
plt.grid(ls=":")
plt.title("$criterion: n=$n")
end
end
plt.tight_layout()
end
plotAllTests (generic function with 1 method)
plotAllTests("WAIC")
plotAllTests("WBIC")
plotAllTests("FreeEnergy")
@show count((WAIC_normal1 .< WAIC_mixnormal) .& (GL_normal1 .< GL_mixnormal))
@show count((WAIC_normal1 .< WAIC_mixnormal) .& (GL_normal1 .≥ GL_mixnormal))
@show count((WAIC_normal1 .> WAIC_mixnormal) .& (GL_normal1 .≤ GL_mixnormal))
@show count((WAIC_normal1 .> WAIC_mixnormal) .& (GL_normal1 .> GL_mixnormal))
;
count((WAIC_normal1 .< WAIC_mixnormal) .& (GL_normal1 .< GL_mixnormal)) = 630 count((WAIC_normal1 .< WAIC_mixnormal) .& (GL_normal1 .≥ GL_mixnormal)) = 182 count((WAIC_normal1 .> WAIC_mixnormal) .& (GL_normal1 .≤ GL_mixnormal)) = 175 count((WAIC_normal1 .> WAIC_mixnormal) .& (GL_normal1 .> GL_mixnormal)) = 13
function show2x2(criterion1::String, criterion2::String, model1::String, model2::String)
C11 = s(criterion1, model1)
C12 = s(criterion1, model2)
C21 = s(criterion2, model1)
C22 = s(criterion2, model2)
count11 = count((C11 .< C12) .& (C21 .< C22))
count12 = count((C11 .< C12) .& (C21 .≥ C22))
count21 = count((C11 .> C12) .& (C21 .≤ C22))
count22 = count((C11 .> C12) .& (C21 .> C22))
@printf("%10s %10s|%-10s selects |\n", "", "", criterion2)
@printf("%10s %10s|----------+----------+\n", "", "")
@printf("%10s %10s|%10s|%10s|\n", "", "", model1, model2)
println("----------+----------+----------+----------+----------")
@printf("%-10s|%10s| %5d | %5d | %5d\n", criterion1, model1, count11, count12, count11+count12)
@printf( "%10s|%10s| %5d | %5d | %5d\n", " selects", model2, count21, count22, count21+count22)
println("----------+----------+----------+----------+----------")
@printf("%10s %10s| %5d | %5d | %5d\n", "", "", count11+count21, count12+count22, length(C11))
end
function show2x2all(criterion1::String, criterion2::String)
model = ["normal1", "mixnormal", "normal"]
println("==========="^5)
for i in 1:3
for j in i+1:3
show2x2(criterion1, criterion2, model[i], model[j])
println("==========="^5)
end
end
end
show2x2all (generic function with 1 method)
show2x2all("WAIC", "GL")
======================================================= |GL selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- WAIC | normal1| 630 | 182 | 812 selects| mixnormal| 175 | 13 | 188 ----------+----------+----------+----------+---------- | 805 | 195 | 1000 ======================================================= |GL selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- WAIC | normal1| 808 | 36 | 844 selects| normal| 156 | 0 | 156 ----------+----------+----------+----------+---------- | 964 | 36 | 1000 ======================================================= |GL selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- WAIC | mixnormal| 658 | 173 | 831 selects| normal| 169 | 0 | 169 ----------+----------+----------+----------+---------- | 827 | 173 | 1000 =======================================================
show2x2all("LOOCV", "GL")
======================================================= |GL selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- LOOCV | normal1| 629 | 183 | 812 selects| mixnormal| 176 | 12 | 188 ----------+----------+----------+----------+---------- | 805 | 195 | 1000 ======================================================= |GL selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- LOOCV | normal1| 809 | 36 | 845 selects| normal| 155 | 0 | 155 ----------+----------+----------+----------+---------- | 964 | 36 | 1000 ======================================================= |GL selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- LOOCV | mixnormal| 666 | 173 | 839 selects| normal| 161 | 0 | 161 ----------+----------+----------+----------+---------- | 827 | 173 | 1000 =======================================================
show2x2all("LOOCV", "WAIC")
======================================================= |WAIC selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- LOOCV | normal1| 810 | 2 | 812 selects| mixnormal| 2 | 186 | 188 ----------+----------+----------+----------+---------- | 812 | 188 | 1000 ======================================================= |WAIC selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- LOOCV | normal1| 844 | 1 | 845 selects| normal| 0 | 155 | 155 ----------+----------+----------+----------+---------- | 844 | 156 | 1000 ======================================================= |WAIC selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- LOOCV | mixnormal| 831 | 8 | 839 selects| normal| 0 | 161 | 161 ----------+----------+----------+----------+---------- | 831 | 169 | 1000 =======================================================
show2x2all("WBIC", "FreeEnergy")
======================================================= |FreeEnergy selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- WBIC | normal1| 693 | 12 | 705 selects| mixnormal| 19 | 276 | 295 ----------+----------+----------+----------+---------- | 712 | 288 | 1000 ======================================================= |FreeEnergy selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- WBIC | normal1| 839 | 28 | 867 selects| normal| 3 | 130 | 133 ----------+----------+----------+----------+---------- | 842 | 158 | 1000 ======================================================= |FreeEnergy selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- WBIC | mixnormal| 780 | 34 | 814 selects| normal| 10 | 176 | 186 ----------+----------+----------+----------+---------- | 790 | 210 | 1000 =======================================================
show2x2all("WBIC", "GL")
======================================================= |GL selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- WBIC | normal1| 549 | 156 | 705 selects| mixnormal| 256 | 39 | 295 ----------+----------+----------+----------+---------- | 805 | 195 | 1000 ======================================================= |GL selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- WBIC | normal1| 831 | 36 | 867 selects| normal| 133 | 0 | 133 ----------+----------+----------+----------+---------- | 964 | 36 | 1000 ======================================================= |GL selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- WBIC | mixnormal| 654 | 160 | 814 selects| normal| 173 | 13 | 186 ----------+----------+----------+----------+---------- | 827 | 173 | 1000 =======================================================
show2x2all("WBIC", "WAIC")
======================================================= |WAIC selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- WBIC | normal1| 660 | 45 | 705 selects| mixnormal| 152 | 143 | 295 ----------+----------+----------+----------+---------- | 812 | 188 | 1000 ======================================================= |WAIC selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- WBIC | normal1| 798 | 69 | 867 selects| normal| 46 | 87 | 133 ----------+----------+----------+----------+---------- | 844 | 156 | 1000 ======================================================= |WAIC selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- WBIC | mixnormal| 760 | 54 | 814 selects| normal| 71 | 115 | 186 ----------+----------+----------+----------+---------- | 831 | 169 | 1000 =======================================================
show2x2all("FreeEnergy", "GL")
======================================================= |GL selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- FreeEnergy| normal1| 552 | 160 | 712 selects| mixnormal| 253 | 35 | 288 ----------+----------+----------+----------+---------- | 805 | 195 | 1000 ======================================================= |GL selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- FreeEnergy| normal1| 806 | 36 | 842 selects| normal| 158 | 0 | 158 ----------+----------+----------+----------+---------- | 964 | 36 | 1000 ======================================================= |GL selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- FreeEnergy| mixnormal| 627 | 163 | 790 selects| normal| 200 | 10 | 210 ----------+----------+----------+----------+---------- | 827 | 173 | 1000 =======================================================
show2x2all("FreeEnergy", "WAIC")
======================================================= |WAIC selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- FreeEnergy| normal1| 674 | 38 | 712 selects| mixnormal| 138 | 150 | 288 ----------+----------+----------+----------+---------- | 812 | 188 | 1000 ======================================================= |WAIC selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- FreeEnergy| normal1| 793 | 49 | 842 selects| normal| 51 | 107 | 158 ----------+----------+----------+----------+---------- | 844 | 156 | 1000 ======================================================= |WAIC selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- FreeEnergy| mixnormal| 754 | 36 | 790 selects| normal| 77 | 133 | 210 ----------+----------+----------+----------+---------- | 831 | 169 | 1000 =======================================================
show2x2all("WAIC", "GL")
======================================================= |GL selects | |----------+----------+ | normal1| mixnormal| ----------+----------+----------+----------+---------- WAIC | normal1| 630 | 182 | 812 selects| mixnormal| 175 | 13 | 188 ----------+----------+----------+----------+---------- | 805 | 195 | 1000 ======================================================= |GL selects | |----------+----------+ | normal1| normal| ----------+----------+----------+----------+---------- WAIC | normal1| 808 | 36 | 844 selects| normal| 156 | 0 | 156 ----------+----------+----------+----------+---------- | 964 | 36 | 1000 ======================================================= |GL selects | |----------+----------+ | mixnormal| normal| ----------+----------+----------+----------+---------- WAIC | mixnormal| 658 | 173 | 831 selects| normal| 169 | 0 | 169 ----------+----------+----------+----------+---------- | 827 | 173 | 1000 =======================================================
function plotchain(chain, paramname)
K = size(chain,2)
m = 0
plt.figure(figsize=(8,8))
for j in 1:K
for i in 1:K
m += 1
fig = plt.subplot(K,K,m)
fig[:set_facecolor]("k")
if i == j
xmin = minimum(chain[:,i])
xmax = maximum(chain[:,i])
x = linspace(xmin, xmax, 200)
ik = InterpKDE(kde(chain[:,i]))
plt.plot(x, pdf.(ik,x), color="pink", lw=1)
plt.grid(ls=":")
plt.xticks(fontsize=6)
plt.yticks(fontsize=6)
j == 3 && plt.xlabel("$(paramname[i])", fontsize=10)
i == 1 && plt.ylabel("$(paramname[i])", fontsize=10)
else
ik = InterpKDE(kde((chain[:,i], chain[:,j])))
plt.scatter(chain[:,i], chain[:,j], c=pdf.(ik, chain[:,i], chain[:,j]), s=2, cmap="CMRmap")
plt.grid(ls=":")
plt.xticks(fontsize=6)
plt.yticks(fontsize=6)
j == 3 && plt.xlabel("$(paramname[i])", fontsize=10)
i == 1 && plt.ylabel("$(paramname[j])", fontsize=10)
end
end
end
plt.tight_layout()
end
modelname = "mixnormal"
chain = eval(Symbol(:Chain_, modelname))[7]
paramname = ["a", "b", "c"]
println("Model: $modelname (n=$n)")
sleep(0.1)
plotchain(chain, paramname)
Model: mixnormal (n=32)
@everywhere mixnormal(a,b,c) = MixtureModel(Normal[Normal(b, 1.0), Normal(c, 1.0)], [1.0-a, a])
mixnormal(w::AbstractVector) = mixnormal(w[1], w[2], w[3])
unlink_mixnormal(w) = [logit(w[1]), w[2], w[3]] # unlink_mixnormal : (0,1)×R^2 -> R^3
link_mixnormal(z) = [invlogit(z[1]), z[2], z[3]] # link_mixnormal : R^3 → (0,1)×R^2
@everywhere normal(mu, sigma) = Normal(mu, sigma)
normal(w::AbstractVector) = normal(w[1], w[2])
unlink_normal(w) = [w[1], log(w[2])] # unlink_normal : R×(0,∞) -> R^2
link_normal(z) = [z[2], exp(z[2])] # link_normal : R^2 → R×(0,∞)
@everywhere normal1(mu::Real) = Normal(mu, 1.0)
normal1(w::AbstractVector) = normal1(w[1])
unlink_normal1(w) = w
link_normal1(z) = z
function makefunc_pred_Bayes(dist_model, chain)
L = size(chain, 1)
function pred_Bayes(x)
return @sum(pdf(dist_model(@view chain[l,:]), x), l, 1:L)/L
end
end
function plotpredictives(t;
show_true=true, show_sample=true, show_mixnormal=true, show_normal1=true, show_normal=true)
Y = Sample[t]
show_mixnormal && chain_mixnormal = Chain_mixnormal[t]
show_normal1 && chain_normal1 = Chain_normal1[t]
show_normal && chain_normal = Chain_normal[t]
show_sample && kde_sample = makefunc_pdfkde(Y)
show_mixnormal && pred_mixnormal = makefunc_pred_Bayes(mixnormal, chain_mixnormal)
show_normal1 && pred_normal1 = makefunc_pred_Bayes(normal1, chain_normal1)
show_normal && pred_normal = makefunc_pred_Bayes(normal, chain_normal)
x = linspace(-4.0, 4.0, 401)
show_true && plt.plot(x, pdf.(dist_true, x), label="true dist", ls="-", color="#1f77b4")
show_sample && plt.plot(x, kde_sample.(x), label="sample KDE", ls="-", color="#9467bd")
show_normal1 && plt.plot(x, pred_normal1.(x), label="normal1", ls="--", color="#ff7f0e")
show_mixnormal && plt.plot(x, pred_mixnormal.(x), label="mixnormal", ls="-.", color="#2ca02c")
show_normal && plt.plot(x, pred_normal.(x), label="normal", ls=":", color="#d62728")
plt.xlabel("x")
plt.ylabel("probability density")
plt.grid(ls=":")
plt.legend(fontsize=8)
plt.title("Predictive PDFs($t): n = $(length(Y))")
end
function plotcomparison(t; show_mixnormal=true, show_normal1=true, show_normal=true)
plt.figure(figsize=(10,3.5))
plt.subplot(121)
plotpredictives(t, show_true=false, show_sample=true,
show_mixnormal=show_mixnormal, show_normal1=show_normal1, show_normal=show_normal)
plt.subplot(122)
plotpredictives(t, show_true=true, show_sample=false,
show_mixnormal=show_mixnormal, show_normal1=show_normal1, show_normal=show_normal)
plt.tight_layout()
end
plotcomparison(1, show_normal=false)
@show t11 = [t for t in 1:Nsims if WAIC_normal1[t] < WAIC_mixnormal[t] && GL_normal1[t] < GL_mixnormal[t]];
@show t12 = [t for t in 1:Nsims if WAIC_normal1[t] < WAIC_mixnormal[t] && GL_normal1[t] > GL_mixnormal[t]];
@show t21 = [t for t in 1:Nsims if WAIC_normal1[t] > WAIC_mixnormal[t] && GL_normal1[t] < GL_mixnormal[t]];
@show t22 = [t for t in 1:Nsims if WAIC_normal1[t] > WAIC_mixnormal[t] && GL_normal1[t] > GL_mixnormal[t]];
t11 = [t for t = 1:Nsims if WAIC_normal1[t] < WAIC_mixnormal[t] && GL_normal1[t] < GL_mixnormal[t]] = [1, 3, 4, 5, 6, 8, 9, 10, 11, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 26, 28, 29, 33, 34, 36, 37, 38, 39, 44, 45, 47, 49, 50, 52, 53, 56, 57, 58, 59, 61, 63, 64, 66, 67, 69, 70, 72, 73, 74, 75, 79, 80, 81, 84, 85, 86, 87, 88, 89, 90, 92, 93, 95, 96, 97, 100, 103, 104, 105, 108, 109, 114, 115, 116, 119, 120, 124, 125, 127, 128, 129, 131, 132, 133, 134, 139, 142, 144, 145, 147, 148, 149, 150, 151, 152, 153, 159, 164, 166, 169, 170, 173, 174, 175, 176, 180, 182, 183, 184, 186, 188, 189, 190, 191, 192, 195, 196, 197, 202, 203, 204, 205, 206, 207, 208, 209, 210, 212, 214, 216, 219, 220, 222, 223, 224, 225, 226, 227, 228, 230, 231, 232, 234, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 253, 254, 256, 257, 258, 260, 261, 262, 263, 264, 265, 266, 267, 268, 270, 273, 275, 276, 277, 278, 281, 282, 283, 284, 285, 286, 287, 288, 290, 291, 293, 294, 295, 296, 297, 299, 301, 304, 306, 311, 312, 314, 315, 316, 318, 319, 320, 321, 323, 327, 328, 329, 330, 333, 334, 335, 338, 339, 340, 341, 343, 345, 346, 349, 352, 354, 355, 356, 357, 359, 360, 365, 369, 371, 372, 374, 376, 377, 379, 380, 383, 384, 387, 388, 390, 391, 392, 393, 394, 395, 397, 399, 402, 403, 404, 405, 407, 409, 410, 416, 417, 418, 419, 420, 421, 423, 425, 430, 431, 432, 434, 435, 436, 438, 439, 442, 443, 444, 445, 446, 447, 450, 452, 456, 457, 458, 460, 461, 463, 464, 465, 468, 469, 470, 472, 475, 476, 477, 478, 480, 483, 484, 485, 486, 488, 490, 492, 493, 497, 498, 499, 500, 501, 502, 503, 505, 506, 509, 510, 511, 512, 513, 514, 515, 516, 517, 519, 521, 523, 524, 525, 529, 532, 533, 534, 535, 536, 538, 539, 540, 541, 543, 544, 547, 551, 552, 553, 554, 556, 557, 558, 559, 561, 562, 565, 567, 568, 569, 570, 572, 573, 574, 575, 576, 578, 580, 582, 584, 585, 586, 588, 593, 594, 596, 597, 598, 599, 601, 602, 603, 605, 606, 607, 608, 610, 611, 613, 614, 615, 617, 618, 622, 625, 630, 631, 632, 633, 634, 636, 637, 638, 639, 640, 641, 644, 647, 648, 649, 653, 654, 655, 656, 659, 660, 661, 663, 664, 666, 668, 669, 670, 672, 673, 675, 676, 679, 680, 681, 682, 683, 686, 688, 690, 692, 695, 696, 697, 700, 701, 702, 704, 707, 712, 713, 714, 716, 717, 718, 721, 722, 724, 727, 730, 732, 733, 734, 735, 738, 739, 742, 745, 750, 751, 754, 755, 756, 757, 759, 761, 764, 765, 767, 769, 770, 771, 772, 775, 776, 777, 778, 779, 780, 781, 783, 784, 785, 786, 789, 791, 792, 793, 795, 797, 798, 800, 802, 803, 806, 807, 808, 809, 814, 817, 818, 819, 820, 821, 823, 824, 828, 829, 830, 832, 833, 834, 835, 837, 838, 839, 841, 842, 843, 844, 845, 849, 850, 851, 852, 856, 858, 860, 863, 864, 870, 871, 873, 874, 875, 876, 877, 879, 881, 882, 883, 884, 885, 886, 888, 889, 890, 891, 892, 893, 894, 897, 898, 900, 901, 902, 904, 906, 907, 908, 909, 910, 911, 913, 914, 915, 916, 918, 920, 922, 923, 924, 925, 927, 930, 931, 932, 933, 936, 937, 938, 939, 940, 941, 942, 943, 944, 946, 948, 953, 954, 956, 957, 958, 959, 962, 963, 964, 965, 966, 967, 968, 969, 970, 972, 973, 974, 975, 976, 977, 979, 980, 982, 983, 985, 986, 987, 988, 990, 991, 992, 996, 997, 998, 999, 1000] t12 = [t for t = 1:Nsims if WAIC_normal1[t] < WAIC_mixnormal[t] && GL_normal1[t] > GL_mixnormal[t]] = [2, 20, 27, 40, 41, 42, 43, 48, 51, 76, 82, 91, 94, 98, 99, 102, 107, 113, 117, 118, 121, 122, 123, 130, 135, 137, 143, 154, 157, 160, 161, 163, 165, 168, 171, 177, 178, 179, 181, 185, 187, 198, 200, 213, 215, 221, 229, 235, 252, 255, 269, 271, 289, 292, 303, 305, 310, 317, 324, 326, 337, 348, 350, 353, 362, 363, 364, 367, 375, 378, 381, 385, 406, 408, 413, 422, 424, 426, 427, 429, 440, 448, 449, 451, 453, 459, 473, 479, 482, 487, 491, 496, 504, 507, 526, 527, 531, 545, 548, 549, 550, 560, 563, 564, 566, 571, 579, 591, 592, 604, 620, 623, 626, 627, 645, 646, 650, 665, 678, 684, 685, 687, 689, 693, 694, 699, 705, 706, 708, 709, 725, 729, 731, 736, 737, 740, 743, 752, 753, 760, 787, 788, 796, 801, 804, 811, 812, 825, 827, 831, 840, 853, 854, 855, 857, 861, 865, 868, 869, 872, 880, 887, 895, 899, 903, 905, 912, 917, 919, 921, 928, 929, 934, 935, 951, 952, 955, 971, 984, 989, 993, 995] t21 = [t for t = 1:Nsims if WAIC_normal1[t] > WAIC_mixnormal[t] && GL_normal1[t] < GL_mixnormal[t]] = [7, 12, 30, 31, 32, 35, 46, 54, 55, 60, 62, 65, 68, 71, 77, 78, 83, 101, 106, 110, 111, 112, 126, 136, 138, 140, 141, 146, 155, 156, 158, 162, 167, 172, 193, 194, 199, 201, 217, 218, 233, 259, 272, 274, 279, 280, 298, 300, 302, 307, 309, 313, 322, 325, 331, 332, 336, 342, 344, 347, 351, 358, 361, 366, 368, 370, 373, 382, 386, 389, 396, 398, 400, 401, 411, 414, 415, 428, 437, 441, 454, 455, 462, 466, 467, 471, 481, 489, 494, 495, 508, 518, 520, 522, 528, 530, 537, 546, 555, 577, 581, 583, 587, 589, 590, 595, 600, 609, 612, 619, 621, 624, 628, 629, 635, 642, 643, 651, 652, 657, 658, 662, 671, 674, 677, 691, 698, 703, 710, 711, 715, 719, 723, 726, 728, 741, 744, 746, 747, 748, 749, 758, 762, 763, 766, 768, 773, 774, 782, 790, 794, 799, 805, 810, 813, 816, 826, 846, 847, 848, 859, 862, 866, 878, 896, 926, 945, 947, 949, 950, 960, 961, 978, 981, 994] t22 = [t for t = 1:Nsims if WAIC_normal1[t] > WAIC_mixnormal[t] && GL_normal1[t] > GL_mixnormal[t]] = [211, 308, 412, 433, 474, 542, 616, 667, 720, 815, 822, 836, 867]
for m in ["normal1", "mixnormal", "normal"]
eval(parse("KL_$m = GL_$m - Shannon"))
eval(parse("TT_$m = T_$m - T_true"))
eval(parse("WT_$m = WAIC_$m - T_true"))
eval(parse("LT_$m = LOOCV_$m - T_true"))
eval(parse("WBT_$m = WBIC_$m - T_true"))
eval(parse("FET_$m = FreeEnergy_$m - T_true"))
eval(parse("WBS_$m = WBIC_$m - Shannon"))
eval(parse("FES_$m = FreeEnergy_$m - Shannon"))
eval(parse("KW_$m = KL_$m + WT_$m"))
eval(parse("VV_$m = KL_$m - TT_$m"))
eval(parse("lambda_$m = mean(KW_$m)/4"))
eval(parse("nu_$m = mean(V_$m)/4"))
eval(parse("nunu_$m = mean(VV_$m)/4"))
eval(parse("KL0_$m = KL_$m - 2lambda_$m"))
eval(parse("WT0_$m = WT_$m - 2lambda_$m"))
end
@show Shannon
println()
@show std(V_normal1), std(V_mixnormal), std(V_normal)
@show 2nu_normal1, 2nu_mixnormal, 2nu_normal
println()
@show std(VV_normal1), std(VV_mixnormal), std(VV_normal)
@show 2nunu_normal1, 2nunu_mixnormal, 2nunu_normal
println()
@show std(KW_normal1), std(KW_mixnormal), std(KW_normal)
@show 2lambda_normal1, 2lambda_mixnormal, 2lambda_normal
;
Shannon = 90.81206612509905 (std(V_normal1), std(V_mixnormal), std(V_normal)) = (0.47302161710052976, 0.8598044435048426, 0.6620573895494153) (2nu_normal1, 2nu_mixnormal, 2nu_normal) = (0.9482629541999565, 1.0742079398050077, 1.80091166751125) (std(VV_normal1), std(VV_mixnormal), std(VV_normal)) = (2.8136299920795556, 2.975801462280639, 3.9827896657255812) (2nunu_normal1, 2nunu_mixnormal, 2nunu_normal) = (0.9852978917341639, 1.0801376952939368, 1.8066087207863344) (std(KW_normal1), std(KW_mixnormal), std(KW_normal)) = (0.30153528059554524, 0.3627806210107444, 0.8367666743839437) (2lambda_normal1, 2lambda_mixnormal, 2lambda_normal) = (0.9330719597667703, 1.155175814734873, 1.9557216808208768)
println("WAIC and GL select normal1.")
println()
@show mean(KL_normal1[t11]), mean(KL_mixnormal[t11])
@show mean(WT_normal1[t11]), mean(WT_mixnormal[t11])
@show mean(KL0_normal1[t11]), mean(KL0_mixnormal[t11])
@show mean(WT0_normal1[t11]), mean(WT0_mixnormal[t11])
sleep(0.1)
@time for i in 1:10
plotcomparison(t11[i], show_normal=false)
plotchain(Chain_mixnormal[t11[i]], ["a","b","c"])
end
WAIC and GL select normal1. (mean(KL_normal1[t11]), mean(KL_mixnormal[t11])) = (0.4496989279611369, 0.6192697549039802) (mean(WT_normal1[t11]), mean(WT_mixnormal[t11])) = (1.3909481805052335, 1.8248345837873012) (mean(KL0_normal1[t11]), mean(KL0_mixnormal[t11])) = (-0.48337303180563274, -0.5359060598308937) (mean(WT0_normal1[t11]), mean(WT0_mixnormal[t11])) = (0.4578762207384633, 0.6696587690524274)
24.064697 seconds (424.22 M allocations: 10.033 GiB, 7.82% gc time)
println("WAIC selects normal1, but GL selects mixnormal.")
println()
@show mean(KL_normal1[t12]), mean(KL_mixnormal[t12])
@show mean(WT_normal1[t12]), mean(WT_mixnormal[t12])
@show mean(KL0_normal1[t12]), mean(KL0_mixnormal[t12])
@show mean(WT0_normal1[t12]), mean(WT0_mixnormal[t12])
sleep(0.1)
@time sleep(0.1)
for i in 1:10
plotcomparison(t12[i], show_normal=false)
plotchain(Chain_mixnormal[t12[i]], ["a","b","c"])
end
WAIC selects normal1, but GL selects mixnormal. (mean(KL_normal1[t12]), mean(KL_mixnormal[t12])) = (2.8122717528102967, 2.6008826627515917) (mean(WT_normal1[t12]), mean(WT_mixnormal[t12])) = (-1.2479196783356956, -0.7378069659802113) (mean(KL0_normal1[t12]), mean(KL0_mixnormal[t12])) = (1.8791997930435276, 1.4457068480167181) (mean(WT0_normal1[t12]), mean(WT0_mixnormal[t12])) = (-2.1809916381024643, -1.8929827807150852)
0.100620 seconds (23 allocations: 1.250 KiB)
println("WAIC selects mixnormal, but GL selects normal1.")
println()
@show mean(KL_normal1[t21]), mean(KL_mixnormal[t21])
@show mean(WT_normal1[t21]), mean(WT_mixnormal[t21])
@show mean(KL0_normal1[t21]), mean(KL0_mixnormal[t21])
@show mean(WT0_normal1[t21]), mean(WT0_mixnormal[t21])
sleep(0.1)
@time for i in 1:10
plotcomparison(t21[i], show_normal=false)
plotchain(Chain_mixnormal[t21[i]], ["a","b","c"])
end
WAIC selects mixnormal, but GL selects normal1. (mean(KL_normal1[t21]), mean(KL_mixnormal[t21])) = (0.7478483900584826, 1.465158638247348) (mean(WT_normal1[t21]), mean(WT_mixnormal[t21])) = (1.5089714082160448, 0.8771456299167988) (mean(KL0_normal1[t21]), mean(KL0_mixnormal[t21])) = (-0.18522356970828766, 0.30998282351247547) (mean(WT0_normal1[t21]), mean(WT0_mixnormal[t21])) = (0.5758994484492743, -0.2780301848180743)
25.108953 seconds (424.22 M allocations: 10.033 GiB, 7.93% gc time)
println("WAIC and GL select mixnormal.")
println()
@show mean(KL_normal1[t22]), mean(KL_mixnormal[t22])
@show mean(WT_normal1[t22]), mean(WT_mixnormal[t22])
@show mean(KL0_normal1[t22]), mean(KL0_mixnormal[t22])
@show mean(WT0_normal1[t22]), mean(WT0_mixnormal[t22])
sleep(0.1)
@time for i in 1:min(length(t22),10)
plotcomparison(t22[i], show_normal=false)
plotchain(Chain_mixnormal[t22[i]], ["a","b","c"])
end
WAIC and GL select mixnormal. (mean(KL_normal1[t22]), mean(KL_mixnormal[t22])) = (3.391511185673285, 3.1694014092475973) (mean(WT_normal1[t22]), mean(WT_mixnormal[t22])) = (-1.323765112803478, -1.5091804590074782) (mean(KL0_normal1[t22]), mean(KL0_mixnormal[t22])) = (2.458439225906514, 2.0142255945127236) (mean(WT0_normal1[t22]), mean(WT0_mixnormal[t22])) = (-2.2568370725702485, -2.6643562737423516)
27.914655 seconds (424.22 M allocations: 10.033 GiB, 7.76% gc time)
@show mean(KL0_normal1), mean(KL0_mixnormal)
@show mean(WT0_normal1), mean(WT0_mixnormal)
println()
println("WAIC and GL select normal1.")
@show mean(KL0_normal1[t11]), mean(KL0_mixnormal[t11])
@show mean(WT0_normal1[t11]), mean(WT0_mixnormal[t11])
@show mean(KL_normal1[t11]), mean(KL_mixnormal[t11])
@show mean(WT_normal1[t11]), mean(WT_mixnormal[t11])
println()
println("WAIC selects normal1, but GL selects mixnormal.")
@show mean(KL0_normal1[t12]), mean(KL0_mixnormal[t12])
@show mean(WT0_normal1[t12]), mean(WT0_mixnormal[t12])
@show mean(KL_normal1[t12]), mean(KL_mixnormal[t12])
@show mean(WT_normal1[t12]), mean(WT_mixnormal[t12])
println()
println("WAIC selects mixnormal, but GL selects normal1.")
@show mean(KL0_normal1[t21]), mean(KL0_mixnormal[t21])
@show mean(WT0_normal1[t21]), mean(WT0_mixnormal[t21])
@show mean(KL_normal1[t21]), mean(KL_mixnormal[t21])
@show mean(WT_normal1[t21]), mean(WT_mixnormal[t21])
println()
println("WAIC and GL select mixnormal.")
@show mean(KL0_normal1[t22]), mean(KL0_mixnormal[t22])
@show mean(WT0_normal1[t22]), mean(WT0_mixnormal[t22])
@show mean(KL_normal1[t22]), mean(KL_mixnormal[t22])
@show mean(WT_normal1[t22]), mean(WT_mixnormal[t22])
;
(mean(KL0_normal1), mean(KL0_mixnormal)) = (0.037034937534206905, 0.005929755488929171) (mean(WT0_normal1), mean(WT0_mixnormal)) = (-0.037034937534207606, -0.005929755488928858) WAIC and GL select normal1. (mean(KL0_normal1[t11]), mean(KL0_mixnormal[t11])) = (-0.48337303180563274, -0.5359060598308937) (mean(WT0_normal1[t11]), mean(WT0_mixnormal[t11])) = (0.4578762207384633, 0.6696587690524274) (mean(KL_normal1[t11]), mean(KL_mixnormal[t11])) = (0.4496989279611369, 0.6192697549039802) (mean(WT_normal1[t11]), mean(WT_mixnormal[t11])) = (1.3909481805052335, 1.8248345837873012) WAIC selects normal1, but GL selects mixnormal. (mean(KL0_normal1[t12]), mean(KL0_mixnormal[t12])) = (1.8791997930435276, 1.4457068480167181) (mean(WT0_normal1[t12]), mean(WT0_mixnormal[t12])) = (-2.1809916381024643, -1.8929827807150852) (mean(KL_normal1[t12]), mean(KL_mixnormal[t12])) = (2.8122717528102967, 2.6008826627515917) (mean(WT_normal1[t12]), mean(WT_mixnormal[t12])) = (-1.2479196783356956, -0.7378069659802113) WAIC selects mixnormal, but GL selects normal1. (mean(KL0_normal1[t21]), mean(KL0_mixnormal[t21])) = (-0.18522356970828766, 0.30998282351247547) (mean(WT0_normal1[t21]), mean(WT0_mixnormal[t21])) = (0.5758994484492743, -0.2780301848180743) (mean(KL_normal1[t21]), mean(KL_mixnormal[t21])) = (0.7478483900584826, 1.465158638247348) (mean(WT_normal1[t21]), mean(WT_mixnormal[t21])) = (1.5089714082160448, 0.8771456299167988) WAIC and GL select mixnormal. (mean(KL0_normal1[t22]), mean(KL0_mixnormal[t22])) = (2.458439225906514, 2.0142255945127236) (mean(WT0_normal1[t22]), mean(WT0_mixnormal[t22])) = (-2.2568370725702485, -2.6643562737423516) (mean(KL_normal1[t22]), mean(KL_mixnormal[t22])) = (3.391511185673285, 3.1694014092475973) (mean(WT_normal1[t22]), mean(WT_mixnormal[t22])) = (-1.323765112803478, -1.5091804590074782)
function plotscatterall(criterion1::String, criterion2::String)
plt.figure(figsize=(10, 3.3))
m = 0
for modelname in ["normal1", "mixnormal", "normal"]
m += 1
fig = plt.subplot(1, 3, m)
plotscatter(criterion1, criterion2, modelname, fig)
end
plt.tight_layout()
end
function plotscatter(criterion1::String, criterion2::String, modelname::String, fig)
C1 = eval(Symbol(criterion1, :_, modelname))
C2 = eval(Symbol(criterion2, :_, modelname))
f = makefunc_pdfkde(C1,C2)
fig[:set_facecolor]("darkgray")
plt.scatter(C1, C2, c=f.(C1,C2), s=3, cmap="CMRmap", alpha=0.5)
plt.xlabel(criterion1)
plt.ylabel(criterion2)
plt.title("$modelname (n = $n)")
x = [minimum(C1), maximum(C1)]
plt.plot(x, x, ls="--", color="k", lw=0.6)
plt.grid(ls=":", color="black", lw=0.3)
end
plotscatterall("WAIC", "LOOCV")
plotscatterall("WT", "LT")
plotscatterall("GL", "WAIC")
plotscatterall("KL", "WT")
plotscatterall("GL", "LOOCV")
plotscatterall("KL", "LT")
plotscatterall("GL", "FreeEnergy")
plotscatterall("KL", "FET")
plotscatterall("GL", "WBIC")
plotscatterall("KL", "WBT")
plotscatterall("FreeEnergy", "WAIC")
plotscatterall("FreeEnergy", "WBIC")
plotscatterall("FET", "WBT")
;
function plotdiffscatterall(criterion1::String, criterion2::String)
plt.figure(figsize=(10, 3.3))
m = 0
for modelname in ["normal1", "mixnormal", "normal"]
m += 1
fig = plt.subplot(1, 3, m)
plotdiffscatter(criterion1, criterion2, modelname, fig)
end
plt.tight_layout()
end
function plotdiffscatter(criterion1::String, criterion2::String, modelname::String, fig)
C1 = eval(Symbol(criterion1, :_, modelname))
C2 = eval(Symbol(criterion2, :_, modelname))
f = makefunc_pdfkde(C1,C2)
fig[:set_facecolor]("darkgray")
plt.scatter(C1, C2-C1, c=f.(C1,C2), s=3, cmap="CMRmap", alpha=0.5)
plt.xlabel(criterion1)
plt.ylabel("$criterion2 \$-\$ $criterion1")
plt.title("$modelname (n = $n)")
x = [minimum(C1), maximum(C1)]
plt.plot(x, x-x, ls="--", color="k", lw=0.6)
plt.grid(ls=":", color="black", lw=0.3)
end
plotdiffscatterall("WAIC", "LOOCV")
plotdiffscatterall("WT", "LT")
plotdiffscatterall("GL", "WAIC")
plotdiffscatterall("KL", "WT")
plotdiffscatterall("GL", "LOOCV")
plotdiffscatterall("KL", "LT")
plotdiffscatterall("GL", "FreeEnergy")
plotdiffscatterall("KL", "FET")
plotdiffscatterall("GL", "WBIC")
plotdiffscatterall("KL", "WBT")
plotdiffscatterall("FreeEnergy", "WAIC")
plotdiffscatterall("FreeEnergy", "WBIC")
plotdiffscatterall("FET", "WBT")
;
println("========== n = $n")
println("mean (std) of T_true - Shannon = $(mean(T_true - Shannon)) ($(std(T_true - Shannon)))")
println("mean (std) of (T_true - Shannon)/(2n) = $(mean(T_true - Shannon)/(2n)) ($(std(T_true - Shannon)/(2n)))")
println()
for modelname in ["normal1", "mixnormal", "normal"]
WT = eval(Symbol(:WT, :_, modelname))
WG = eval(Symbol(:WAIC, :_, modelname)) - eval(Symbol(:GL, :_, modelname))
LT = eval(Symbol(:LT, :_, modelname))
LG = eval(Symbol(:LOOCV, :_, modelname)) - eval(Symbol(:GL, :_, modelname))
FET = eval(Symbol(:FET, :_, modelname))
FES = eval(Symbol(:FES, :_, modelname))
WBT = eval(Symbol(:WBT, :_, modelname))
WBS = eval(Symbol(:WBS, :_, modelname))
println("========== $modelname (n = $n)")
println("mean (std) of WAIC - T_true = $(mean(WT)) ($(std(WT)))")
println("mean (std) of LOOCV - T_true = $(mean(LT)) ($(std(LT)))")
println("mean (std) of FreeEnergy - T_true = $(mean(FET)) ($(std(FET)))")
println("mean (std) of WBIC - T_true = $(mean(WBT)) ($(std(WBT)))")
println("mean (std) of WAIC - GL = $(mean(WG)) ($(std(WG)))")
println("mean (std) of LOOCV - GL = $(mean(LG)) ($(std(LG)))")
println("mean (std) of FreeEnergy - Shannon = $(mean(FES)) ($(std(FES)))")
println("mean (std) of WBIC - Shannon = $(mean(WBS)) ($(std(WBS)))")
println()
end
========== n = 32 mean (std) of T_true - Shannon = -0.0983617194064123 (7.694310808571493) mean (std) of (T_true - Shannon)/(2n) = -0.0015369018657251923 (0.12022360638392958) ========== normal1 (n = 32) mean (std) of WAIC - T_true = 0.8960370222325631 (1.5004116791513726) mean (std) of LOOCV - T_true = 0.8977237017611316 (1.5006163913469779) mean (std) of FreeEnergy - T_true = 1.2609557301661718 (1.7162923302533113) mean (std) of WBIC - T_true = 1.8569652408466168 (2.0681075319960573) mean (std) of WAIC - GL = -0.17243159447482673 (7.9773306007179405) mean (std) of LOOCV - GL = -0.17074491494625815 (7.977763591516943) mean (std) of FreeEnergy - Shannon = 1.1625940107597597 (7.666203176930304) mean (std) of WBIC - Shannon = 1.758603521440205 (7.7665404714558175) ========== mixnormal (n = 32) mean (std) of WAIC - T_true = 1.1492460592459444 (1.5218424853234411) mean (std) of LOOCV - T_true = 1.147950448687549 (1.5235623355264354) mean (std) of FreeEnergy - T_true = 1.692059198221755 (1.7877247594393728) mean (std) of WBIC - T_true = 2.411311971082669 (1.991077828522291) mean (std) of WAIC - GL = -0.11022123038427012 (7.245313316555237) mean (std) of LOOCV - GL = -0.11151684094266495 (7.239324044442848) mean (std) of FreeEnergy - Shannon = 1.5936974788153424 (6.962418987983378) mean (std) of WBIC - Shannon = 2.3129502516762566 (6.984221756458584) ========== normal (n = 32) mean (std) of WAIC - T_true = 1.950024627545792 (2.1290880841650597) mean (std) of LOOCV - T_true = 2.0028743741341675 (2.144785562387807) mean (std) of FreeEnergy - T_true = 2.7944988763040546 (2.4243462877296116) mean (std) of WBIC - T_true = 4.202293030525466 (3.052600632531801) mean (std) of WAIC - GL = -0.10975582595658047 (8.558186765597165) mean (std) of LOOCV - GL = -0.05690607936820591 (8.564761570196836) mean (std) of FreeEnergy - Shannon = 2.696137156897642 (8.10096824594339) mean (std) of WBIC - Shannon = 4.103931311119053 (8.245736767533247)
function plotallcomparisonsby(criterion::String)
plt.figure(figsize=(10, 3.3))
fig = plt.subplot(131)
plotcomparisonby(criterion, "normal1", "mixnormal", fig)
fig = plt.subplot(132)
plotcomparisonby(criterion, "mixnormal", "normal", fig)
fig = plt.subplot(133)
plotcomparisonby(criterion, "normal1", "normal", fig)
plt.tight_layout()
end
function plotcomparisonby(criterion::String, modelname1::String, modelname2::String, fig)
C1 = eval(Symbol(criterion, :_, modelname1))
C2 = eval(Symbol(criterion, :_, modelname2))
f = makefunc_pdfkde(C1,C2)
fig[:set_facecolor]("darkgray")
plt.scatter(C1, C2, c=f.(C1,C2), s=3, cmap="CMRmap", alpha=0.5)
plt.xlabel(modelname1)
plt.ylabel(modelname2)
plt.title("$criterion (n = $n)")
x = [minimum(C1), maximum(C1)]
plt.plot(x, x, ls="--", color="k", lw=0.6)
plt.grid(ls=":", color="black", lw=0.3)
end
#plotallcomparisonsby("GL")
plotallcomparisonsby("KL")
#plotallcomparisonsby("WAIC")
plotallcomparisonsby("WT")
#plotallcomparisonsby("LOOCV")
plotallcomparisonsby("LT")
#plotallcomparisonsby("FreeEnergy")
plotallcomparisonsby("FET")
#plotallcomparisonsby("WBIC")
plotallcomparisonsby("WBT")
;
function plotallcomparisonsby(criterion1::String, criterion2::String)
plt.figure(figsize=(10, 6.6))
fig1 = plt.subplot(231)
fig2 = plt.subplot(232)
fig3 = plt.subplot(233)
fig4 = plt.subplot(234)
fig5 = plt.subplot(235)
fig6 = plt.subplot(236)
plotcomparisonby(criterion1, criterion2, "normal1", "mixnormal", fig1, fig4)
plotcomparisonby(criterion1, criterion2, "mixnormal", "normal", fig2, fig5)
plotcomparisonby(criterion1, criterion2, "normal1", "normal", fig3, fig6)
plt.tight_layout()
end
function plotcomparisonby(
criterion1::String, criterion2::String,
modelname1::String, modelname2::String,
fig1, fig2)
C11 = eval(Symbol(criterion1, :_, modelname1))
C12 = eval(Symbol(criterion1, :_, modelname2))
C21 = eval(Symbol(criterion2, :_, modelname1))
C22 = eval(Symbol(criterion2, :_, modelname2))
C11_11 = C11[(C11 .< C12) .& (C21 .< C22)]
C12_11 = C12[(C11 .< C12) .& (C21 .< C22)]
C11_12 = C11[(C11 .< C12) .& (C21 .> C22)]
C12_12 = C12[(C11 .< C12) .& (C21 .> C22)]
C11_21 = C11[(C11 .> C12) .& (C21 .< C22)]
C12_21 = C12[(C11 .> C12) .& (C21 .< C22)]
C11_22 = C11[(C11 .> C12) .& (C21 .> C22)]
C12_22 = C12[(C11 .> C12) .& (C21 .> C22)]
C21_11 = C21[(C11 .< C12) .& (C21 .< C22)]
C22_11 = C22[(C11 .< C12) .& (C21 .< C22)]
C21_12 = C21[(C11 .< C12) .& (C21 .> C22)]
C22_12 = C22[(C11 .< C12) .& (C21 .> C22)]
C21_21 = C21[(C11 .> C12) .& (C21 .< C22)]
C22_21 = C22[(C11 .> C12) .& (C21 .< C22)]
C21_22 = C21[(C11 .> C12) .& (C21 .> C22)]
C22_22 = C22[(C11 .> C12) .& (C21 .> C22)]
fig1[:set_facecolor]("darkgray")
fig1[:scatter](C11_11, C12_11, s=3, color="red", alpha=0.5)
fig1[:scatter](C11_12, C12_12, s=3, color="yellow", alpha=0.5)
fig1[:scatter](C11_21, C12_21, s=3, color="cyan", alpha=0.5)
fig1[:scatter](C11_22, C12_22, s=3, color="blue", alpha=0.5)
fig1[:set_xlabel](modelname1)
fig1[:set_ylabel](modelname2)
fig1[:set_title]("$criterion1 (n = $n)")
x = [minimum(C11), maximum(C11)]
fig1[:plot](x, x, ls="--", color="white", lw=1)
fig1[:grid](ls=":", color="black", lw=0.3)
fig2[:set_facecolor]("darkgray")
fig2[:scatter](C21_11, C22_11, s=3, color="red", alpha=0.5)
fig2[:scatter](C21_12, C22_12, s=3, color="yellow", alpha=0.5)
fig2[:scatter](C21_21, C22_21, s=3, color="cyan", alpha=0.5)
fig2[:scatter](C21_22, C22_22, s=3, color="blue", alpha=0.5)
fig2[:set_xlabel](modelname1)
fig2[:set_ylabel](modelname2)
fig2[:set_title]("$criterion2 (n = $n)")
x = [minimum(C21), maximum(C21)]
fig2[:plot](x, x, ls="--", color="white", lw=1)
fig2[:grid](ls=":", color="black", lw=0.3)
end
plotallcomparisonsby("KL", "WT")
plotallcomparisonsby("KL", "LT")
plotallcomparisonsby("WT", "LT")
plotallcomparisonsby("KL", "FET")
plotallcomparisonsby("KL", "WBT")
plotallcomparisonsby("FET", "WBT")
;
function make_loglik(modelname::String, i)
Y = Sample[i]
chain = eval(Symbol(:Chain_, modelname))[i]
dist_model = eval(Symbol(modelname))
n = length(Y)
L = size(chain, 1)
loglik = Array{Float64,2}(L,n)
for k in 1:n
for l in 1:L
loglik[l,k] = logpdf(dist_model(chain[l,:]), Y[k])
end
end
return loglik
end
function FreeEnergyErrorE2nLn(loglik)
E2nLn = makefunc_E2nLn(loglik)
FreeEnergy, Error = quadgk(E2nLn, 0.0, 1.0)
return FreeEnergy, Error, E2nLn
end
function makefunc_E2nLn(loglik)
L = size(loglik)[1]
negloglik = -sum(loglik, 2)
negloglik_n = negloglik .- maximum(negloglik)
function E2nLn(beta)
Edenominator = @sum( exp((1-beta)*negloglik_n[l]), l, 1:L)/L
if Edenominator == zero(Edenominator) || !isfinite(Edenominator)
return zero(Edenominator)
end
Enumerator = @sum(negloglik[l]*exp((1-beta)*negloglik_n[l]), l, 1:L)/L
return 2*Enumerator/Edenominator
end
return E2nLn
end
function plotFreeEnergy(modelname::String, i)
loglik = make_loglik(modelname, i)
FreeEnergy, Error, E2nLn = FreeEnergyErrorE2nLn(loglik)
WBIC = eval(Symbol(:WBIC_, modelname))[i]
#@show WBIC, FreeEnergy, Error, WBIC
#sleep(0.1)
#plt.figure(figsize=(4,3))
x = linspace(0,1,101)
plt.plot(x, E2nLn.(x), label="\$2E[nL_n|\\beta]\$", color="blue")
plt.plot([0,1], fill(FreeEnergy,2), label="FreeEnergy", ls="--", color="blue")
plt.plot([0,1], fill(WBIC,2), label="WBIC", ls="-.", color="red")
plt.axvline(1/log(n), ls="--", label="x=1/log(n)", ls=":", color="red")
plt.grid(ls=":")
plt.legend(fontsize=9)
plt.title("$modelname $i (n=$n)")
end
function plotFreeEnergyAll(i)
plt.figure(figsize=(10,3.3))
plt.subplot(131)
plotFreeEnergy("normal", i)
plt.subplot(132)
plotFreeEnergy("mixnormal", i)
plt.subplot(133)
plotFreeEnergy("normal1", i)
plt.tight_layout()
end
@time for i in 1:10
plotFreeEnergyAll(i)
end
5.198741 seconds (34.46 M allocations: 854.868 MiB, 4.16% gc time)
@time for i in 1:10
plotFreeEnergyAll(rand(1:Nsims))
end
4.535613 seconds (34.37 M allocations: 850.453 MiB, 3.95% gc time)
function Z_exactf(beta, x)
n = length(x)
return 1/sqrt( (2*pi)^(n*beta) * (n*beta+1) )*exp(-( beta*sum(x.^2) - beta^2/(n*beta+1)*sum(x)^2 )/2)
end
function F_exactf(beta, x)
n = length(x)
return beta*sum(x.^2) - beta^2/(n*beta+1)*sum(x)^2 + log(n*beta+1) + n*beta*log(2*pi)
end
function E2nLn_exactf(beta, x)
n = length(x)
return sum(x.^2) - ( 2*beta/(n*beta+1) - n*beta^2/(n*beta+1)^2 )*sum(x)^2 + n/(n*beta+1) + n*log(2*pi)
end
function Z_quadgk(beta, x)
n = length(x)
f(mu) = 1/sqrt((2*pi)^(n*beta+1))*exp(-( beta*sum((mu.-x).^2) + mu^2 )/2)
return quadgk(f, -20, 20)[1]
end
function F_quadgk(beta, x)
f(beta) = E2nLn_exactf(beta, x)
return quadgk(f, 0, beta)[1]
end
beta = 0.7
x = rand(n)
@show -2*log(Z_exactf(beta,x))
@show F_exactf(beta,x)
@show -2*log(Z_quadgk(beta,x))
@show F_quadgk(beta,x)
println()
Z_normal1(beta, x) = Z_exactf(beta, x)
F_normal1(beta, x) = F_exactf(beta, x)
E2nLn_normal1(beta, x) = E2nLn_exactf(beta, x)
FreeEnergyExact_normal1 = Array{Float64}(Nsims)
for i in 1:Nsims
FreeEnergyExact_normal1[i] = F_normal1(1.0, Sample[i])
end
@show mean(FreeEnergy_normal1 - FreeEnergyExact_normal1)
@show std(FreeEnergy_normal1 - FreeEnergyExact_normal1)
@show mean(WBIC_normal1 - FreeEnergyExact_normal1)
@show std(WBIC_normal1 - FreeEnergyExact_normal1)
sleep(0.1)
FEET_normal1 = FreeEnergyExact_normal1 - T_true
x = [minimum(FEET_normal1), maximum(FEET_normal1)]
ymin = min(minimum(FET_normal1), minimum(WBT_normal1))
ymax = max(maximum(FET_normal1), maximum(WBT_normal1))
f1 = makefunc_pdfkde(FEET_normal1, FET_normal1)
c1 = f1.(FEET_normal1, FET_normal1)
f2 = makefunc_pdfkde(FEET_normal1, WBT_normal1)
c2 = f2.(FEET_normal1, WBT_normal1)
plt.figure(figsize=(8, 4.2))
fig1 = plt.subplot(121)
fig1[:set_facecolor]("darkgray")
plt.scatter(FEET_normal1, FET_normal1, c=c1, label="FreeEnergy \$-\$ T_true", s=5, alpha=0.5, cmap="CMRmap")
plt.plot(x, x, ls=":", color="white")
plt.xlabel("Exact Free Energy \$-\$ T_true")
plt.ylim(ymin, ymax)
plt.grid(ls=":", color="lightgray")
plt.legend()
plt.title("normal1 (n=$n)")
fig2 = plt.subplot(122)
fig2[:set_facecolor]("darkgray")
plt.scatter(FEET_normal1, WBT_normal1, c=c2, label="WBIC \$-\$ T_true", s=5, alpha=0.5, cmap="CMRmap")
plt.plot(x, x, ls=":", color="white")
plt.xlabel("Exact Free Energy \$-\$ T_true")
plt.ylim(ymin, ymax)
plt.grid(ls=":", color="lightgray")
plt.legend()
plt.title("normal1 (n=$n)")
plt.tight_layout()
;
-2 * log(Z_exactf(beta, x)) = 46.66032071940648 F_exactf(beta, x) = 46.66032071940647 -2 * log(Z_quadgk(beta, x)) = 46.66032071940648 F_quadgk(beta, x) = 46.66032071940647 mean(FreeEnergy_normal1 - FreeEnergyExact_normal1) = -1.2173034602523751 std(FreeEnergy_normal1 - FreeEnergyExact_normal1) = 0.8761541749985949 mean(WBIC_normal1 - FreeEnergyExact_normal1) = -0.6212939495719301 std(WBIC_normal1 - FreeEnergyExact_normal1) = 1.439740113173562
function plotFreeEnergy_normal1(i)
modelname = "normal1"
loglik = make_loglik(modelname, i)
FreeEnergy, Error, E2nLn = FreeEnergyErrorE2nLn(loglik)
WBIC = eval(Symbol(:WBIC_, modelname))[i]
x = linspace(0.01,1,101)
f(beta) = E2nLn_normal1(beta, Sample[i])
plt.plot(x, f.(x), label="exact \$2E[nL_n|\\beta]\$", color="cyan")
plt.plot([0,1], fill(FreeEnergyExact_normal1[i],2), label="FreeEnergyExact", ls="--", color="cyan")
plt.plot(x, E2nLn.(x), label="MCMC \$2E[nL_n|\\beta]\$", color="blue")
plt.plot([0,1], fill(FreeEnergy,2), label="FreeEnergy", ls="--", color="blue")
plt.plot([0,1], fill(WBIC,2), label="WBIC", ls="-.", color="red")
plt.axvline(1/log(n), ls="--", label="x=1/log(n)", ls=":", color="red")
plt.grid(ls=":")
plt.legend(fontsize=9)
plt.title("$modelname $i (n=$n)")
end
i=1
plt.figure(figsize=(10,10))
plt.subplot(221); plotFreeEnergy_normal1(i)
plt.subplot(222); plotFreeEnergy_normal1(i+1)
plt.subplot(223); plotFreeEnergy_normal1(i+1)
plt.subplot(224); plotFreeEnergy_normal1(i+1)
plt.tight_layout()