n = 128 @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 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 #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"); 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("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)) ; 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("WAIC", "GL") show2x2all("LOOCV", "GL") show2x2all("LOOCV", "WAIC") show2x2all("WBIC", "FreeEnergy") show2x2all("WBIC", "GL") show2x2all("WBIC", "WAIC") show2x2all("FreeEnergy", "GL") show2x2all("FreeEnergy", "WAIC") show2x2all("WAIC", "GL") 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) @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]]; 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 ; 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 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 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 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 @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]) ; 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 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 @time for i in 1:10 plotFreeEnergyAll(rand(1:Nsims)) end 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() ; 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="exact FreeEnergy", ls="--", color="cyan") plt.plot(x, E2nLn.(x), label="MCMC \$2E[nL_n|\\beta]\$", color="blue") plt.plot([0,1], fill(FreeEnergy,2), label="MCMC FreeEnergy", ls="--", color="blue") plt.plot([0,1], fill(WBIC,2), label="MCMC 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()