情報量規準のシミュレーション結果のデータの解析 n=128

黒木玄

2017-11-16

このノートブックは

http://nbviewer.jupyter.org/gist/genkuroki/d688c6e5e8c427dba9856d37a12bbee1

で公開されているJupyter notebookで作成したデータファイルを解析するためのノートブックである.

使い方

(1) https://genkuroki.github.io/documents/Jupyter/#2017-11-16 から次の3つのファイルをダウンロードしてこのノートブックと同じ場所に置いておく.

(2) このノートブックを実行する.

諸定義とデータの読み込み

In [1]:
n = 128
Out[1]:
128
In [2]:
@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
 15.064773 seconds (14.84 M allocations: 958.285 MiB, 3.23% gc time)
Out[2]:
@sum (macro with 1 method)
In [3]:
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
Out[3]:
plotAllModelSelections (generic function with 1 method)
In [4]:
#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 = 128)
2λ (std) = 1.3132198826488952 (0.11051169524376306)
2ν (std) = 1.2434679742589019 (0.3893114662850205)
======================================================================
Name                | mean (std)
----------------------------------------------------------------------
T                   | 362.4013518528037 (14.810658604163986)
V                   | 2.4869359485178038 (0.778622932570041)
WAIC = T + V        | 364.8882878013216 (15.517291034838108)
LOOCV               | 364.89090087933636 (15.516582394095808)
GL                  | 364.57523818609155 (1.3700130033861326)
----------------------------------------------------------------------
WAIC - GL           | 0.31304961523019587 (15.217151166016553)
LOOCV - GL          | 0.3156626932450165 (15.216473397507801)
LOOCV - WAIC        | 0.002613078014820303 (0.016397337807532348)
----------------------------------------------------------------------
KL = GL - S         | 1.326973685695199 (1.3700130033861324)
TT = T    - T_true  | -1.1874698689152123 (1.8987780812446018)
WT = WAIC - T_true  | 1.2994660796025912 (1.5205962827214443)
KW = KL + WT        | 2.6264397652977904 (0.22102339048752612)
----------------------------------------------------------------------
WBIC                | 367.24496873183443 (15.095970552245388)
FreeEnergy          | 365.5825765954295 (15.096921251818156)
WBIC - FreeEnergy   | 1.6623921364047805 (1.1796949669415928)
WBIC - T_true       | 3.6561470101153746 (2.7463666443529444)
FreeEnergy - T_true | 1.9937548737105935 (1.9387579307210285)
======================================================================
correlation of TT and KL = -0.9255381674740645
correlation of WT and KL = -0.9937174649007365
======================================================================


========== normal1 model (n = 128)
2λ (std) = 0.9872297188536074 (0.09104830972127553)
2ν (std) = 0.9901284690677659 (0.13409535988613622)
======================================================================
Name                | mean (std)
----------------------------------------------------------------------
T                   | 362.5914240857164 (15.829997087659219)
V                   | 1.9802569381355317 (0.26819071977227243)
WAIC = T + V        | 364.57168102385185 (16.068948210559768)
LOOCV               | 364.5717505913329 (16.06894672700084)
GL                  | 364.2398646359707 (1.363604295771297)
----------------------------------------------------------------------
WAIC - GL           | 0.3318163878813046 (16.18741792385838)
LOOCV - GL          | 0.3318859553623783 (16.18741488768318)
LOOCV - WAIC        | 6.956748107404565e-5 (0.0007035589518467509)
----------------------------------------------------------------------
KL = GL - S         | 0.9916001355743571 (1.3636042957712966)
TT = T    - T_true  | -0.9973976360026738 (1.390112033637219)
WT = WAIC - T_true  | 0.9828593021328578 (1.4012246710907141)
KW = KL + WT        | 1.9744594377072149 (0.18209661944255107)
----------------------------------------------------------------------
WBIC                | 366.141186165668 (16.103590815617803)
FreeEnergy          | 364.95603733863896 (15.999726040888252)
WBIC - FreeEnergy   | 1.1851488270292296 (0.8469179151774902)
WBIC - T_true       | 2.55236444394917 (2.084788968791164)
FreeEnergy - T_true | 1.3672156169199396 (1.558645134078157)
======================================================================
correlation of TT and KL = -0.995186784434587
correlation of WT and KL = -0.9916931960344842
======================================================================


========== normal model (n = 128)
2λ (std) = 1.9912157379734927 (0.2261384408382408)
2ν (std) = 1.934185419939284 (0.20383237518918526)
======================================================================
Name                | mean (std)
----------------------------------------------------------------------
T                   | 361.70670901233996 (16.070556650363827)
V                   | 3.868370839878568 (0.4076647503783705)
WAIC = T + V        | 365.57507985221827 (16.073936379895116)
LOOCV               | 365.5789309425814 (16.07394617064168)
GL                  | 365.24443784584395 (1.9146146999184985)
----------------------------------------------------------------------
WAIC - GL           | 0.3306420063744822 (16.169756406600527)
LOOCV - GL          | 0.33449309673736194 (16.16972509139636)
LOOCV - WAIC        | 0.0038510903628798587 (0.00753071050782785)
----------------------------------------------------------------------
KL = GL - S         | 1.9961733454476547 (1.9146146999184979)
TT = T    - T_true  | -1.8821127093792347 (1.9725632751261004)
WT = WAIC - T_true  | 1.986258130499333 (1.981920621854508)
KW = KL + WT        | 3.9824314759469854 (0.4522768816764816)
----------------------------------------------------------------------
WBIC                | 368.6885413201743 (16.30779423340957)
FreeEnergy          | 366.30698047655505 (16.125020854918912)
WBIC - FreeEnergy   | 2.3815608436194133 (1.4037875138591407)
WBIC - T_true       | 5.099719598455353 (3.260819837669346)
FreeEnergy - T_true | 2.7181587548359407 (2.3081256996668325)
======================================================================
correlation of TT and KL = -0.9793187354190306
correlation of WT and KL = -0.9736436602286695
======================================================================