versioninfo()

@time begin
    using Mamba
    
    using KernelDensity
    function makefunc_pdfkde(X)
        local ik = InterpKDE(kde(X))
        local pdfkde(x) = pdf(ik, x)
        return pdfkde
    end
    function makefunc_pdfkde(X,Y)
        local ik = InterpKDE(kde((X,Y)))
        local 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

using PyCall
@pyimport matplotlib.animation as anim
function showgif(filename)
    open(filename) do f
        base64_video = base64encode(f)
        display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
    end
end

macro sum(f_k, k, itr)
    quote
        begin
            local 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

@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

function sample2model_mixnormal(Y;
        dist_model = mixnormal,
        a0 = 0.5,
        b0 = 0.0,
        c0 = 0.0,
        prior_a = Uniform(0.0, 1.0),
        prior_b = Normal(0.0, 1.0),
        prior_c = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
    )
    
    local model = Model(
        y = Stochastic(1, (a, b, c) -> dist_model(a, b, c), false),
        a = Stochastic(() -> prior_a, true),
        b = Stochastic(() -> prior_b, true),
        c = Stochastic(() -> prior_b, true),
    )
    local scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :a => data[:a0],
            :b => data[:b0],
            :c => data[:c0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

@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,∞)

function sample2model_normal(Y;
        dist_model = normal,
        mu0    = 0.0,
        sigma0 = 1.0,
        prior_mu = Normal(0.0, 1.0),
        prior_sigma = Truncated(Normal(1.0, 1.0), 0, Inf),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :sigma0 => sigma0,
    )
    
    local model = Model(
        y = Stochastic(1, (mu, sigma) -> dist_model(mu, sigma), false),
        mu = Stochastic(() -> prior_mu, true),
        sigma = Stochastic(() -> prior_sigma, true),
    )
    local scheme = [
        NUTS([:mu, :sigma])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
            :sigma => data[:sigma0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

@everywhere normal1(mu::Real) = Normal(mu, 1.0)

normal1(w::AbstractVector) = normal1(w[1])
unlink_normal1(w) = w
link_normal1(z)   = z

function sample2model_normal1(Y;
        dist_model = normal1,
        mu0    = 0.0,
        prior_mu = Normal(0.0, 1.0),
        chains = 2
    )
    local data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
    )
    
    local model = Model(
        y = Stochastic(1, mu -> dist_model(mu), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    local scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    local inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end

## Summary
function showsummary(sim;
        sortkeys=true, figsize_t=(8, 3), figsize_c=(8, 3.5),
        show_describe=true, show_gelman=true, plot_trace=true, plot_contour=true)
    ## Summary of MCMC
    if show_describe
        println("\n========== Summary:\n")
        display(describe(sim))
    end

    # Convergence Diagnostics
    if show_gelman && length(sim.chains) ≥ 2 
       println("========== Gelman Diagnostic:")
       show(gelmandiag(sim))
    end

    ## Plot
    sleep(0.1)
    if plot_trace
        #draw(plot(sim), fmt=:png, width=10inch, height=3.5inch, nrow=1, ncol=2, ask=false)
        pyplot_trace(sim, sortkeys=sortkeys, figsize=figsize_t)
    end
    if plot_contour
        #draw(plot(sim, :contour), fmt=:png, width=10inch, height=4.5inch, nrow=1, ncol=2, ask=false)
        pyplot_contour(sim, sortkeys=sortkeys, figsize=figsize_c)
    end
end

## plot traces
function pyplot_trace(sim; sortkeys = false, figsize = (8, 3))
    if sortkeys
        keys_sim = sort(keys(sim))
    else
        keys_sim = keys(sim)
    end
    for var in keys_sim
        plt.figure(figsize=figsize)
        
        plt.subplot(1,2,1)
        for k in sim.chains
            plt.plot(sim.range, sim[:,var,:].value[:,1,k], label="$k", lw=0.4, alpha=0.8)
        end
        plt.xlabel("iterations")
        plt.grid(ls=":")
        #plt.legend(loc="upper right")
        plt.title("trace of $var", fontsize=10)
        
        plt.subplot(1,2,2)
        local xmin = quantile(vec(sim[:,var,:].value), 0.005)
        local xmax = quantile(vec(sim[:,var,:].value), 0.995)
        for k in sim.chains
            local chain = sim[:,var,:].value[:,1,k]
            local pdfkde = makefunc_pdfkde(chain)
            local x = linspace(xmin, xmax, 201)
            plt.plot(x, pdfkde.(x), label="$k", lw=0.8, alpha=0.8)
        end
        plt.xlabel("$var")
        plt.grid(ls=":")
        plt.title("empirical posterior pdf of $var", fontsize=10)

        plt.tight_layout()
    end
end

# plot contours
function pyplot_contour(sim; sortkeys = true, figsize = (8, 3.5))
    if sortkeys
        keys_sim = sort(keys(sim))
    else
        keys_sim = keys(sim)
    end
    local m = 0
    local K = length(keys_sim)
    for i in 1:K
        for j in i+1:K
            m += 1
            mod1(m,2) == 1 && plt.figure(figsize=figsize)
            plt.subplot(1,2, mod1(m,2))

            local u = keys_sim[i]
            local v = keys_sim[j]
            local X = vec(sim[:,u,:].value)
            local Y = vec(sim[:,v,:].value)
            local pdfkde = makefunc_pdfkde(X,Y)
            local xmin = quantile(X, 0.005)
            local xmax = quantile(X, 0.995)
            local ymin = quantile(Y, 0.005)
            local ymax = quantile(Y, 0.995)
            local x = linspace(xmin, xmax, 200)
            local y = linspace(ymin, ymax, 200)
            
            plt.pcolormesh(x', y, pdfkde.(x',y), cmap="CMRmap")
            plt.colorbar()
            plt.grid(ls=":")
            plt.xlabel(u)
            plt.ylabel(v)
            plt.title("posterior of ($u, $v)", fontsize=10)

            mod1(m,2) == 2 && plt.tight_layout()

            if 2*m == K*(K-1) && mod1(m,2) == 1
                plt.subplot(1,2,2)
                
                plt.pcolormesh(y', x, pdfkde.(x,y'), cmap="CMRmap")
                plt.colorbar()
                plt.grid(ls=":")
                plt.xlabel(v)
                plt.ylabel(u)
                plt.title("posterior of ($v, $u)", fontsize=10)

                plt.tight_layout()
            end
        end
    end
end

# loglik[l,i] = lpdf(w_l, Y_i) と chain[l,:] = w_l を取り出す函数を作る函数
#
function makefunc_loglikchainof(lpdf, symbols, Y)
    #
    # loglikchainof(sim) で loglik[l,i] と chain[l,:] が抽出される
    #
    local function loglikchainof(sim)
        local val = sim[:, symbols, :].value
        local chain = vcat((val[:,:,k] for k in 1:size(val,3))...)
        local L = size(chain,1)
        local n = length(Y)
        local loglik = Array{Float64, 2}(L, n)
        for i in 1:n
            for l in 1:L
                loglik[l,i] = lpdf(chain[l,:], Y[i])
            end
        end
        return loglik, chain
    end
    return loglikchainof
end

# 予測分布函数 p^*(x,y) = mean of { lpdf(w_l, y) }_{l=1}^L を作る函数
#
function makefunc_pdfpred(lpdf, chain)
    local L = size(chain,1)
    local pred_Bayes(y) = @sum(exp(lpdf((@view chain[l,:]), y)), l, 1:L)/L
    return pred_Bayes
end

# loglik[l,i] からWAICを計算する函数
#
function WAICof(loglik)
    local L, n
    L, n = size(loglik)
    local T_n = -mean(log(mean(exp(loglik[l,i]) for l in 1:L)) for i in 1:n)
    local V_n  = sum(var(loglik[:,i], corrected=false) for i in 1:n)
    local WAIC = 2*n*T_n + 2*V_n
    return WAIC, 2*n*T_n, 2*V_n
end

# loglik[l,i] からLOOCVを素朴に計算する函数
#
function LOOCVof(loglik)
    local L, n
    L, n = size(loglik)
    local LOOCV = 2*sum(log(mean(exp(-loglik[l,i]) for l in 1:L)) for i in 1:n)
    return LOOCV
end

# 自由エネルギー(の2倍)を計算するための函数
# 
# 自由エネルギーの2の逆温度微分は E^β_w[2n L_n] に等しいので、
# それを β=0.0 から 1.0 まで数値積分すれば自由エネルギーの2倍を計算できる。
#
function FreeEnergyof(loglik)
    local E2nLn = makefunc_E2nLn(loglik)
    local F = quadgk(E2nLn, 0.0, 1.0)[1]
    return F, E2nLn
end

function makefunc_E2nLn(loglik)
    local L = size(loglik)[1]
    local negloglik = -sum(loglik, 2)
    local negloglik_n = negloglik .- maximum(negloglik)
    local function E2nLn(beta)
        local Edenominator = @sum(             exp((1-beta)*negloglik_n[l]), l, 1:L)/L
        if Edenominator == zero(Edenominator) || !isfinite(Edenominator)
            return zero(Edenominator)
        end
        local Enumerator   = @sum(negloglik[l]*exp((1-beta)*negloglik_n[l]), l, 1:L)/L
        return 2*Enumerator/Edenominator
    end
    return E2nLn
end

# loglik[l,i] からWBICを計算する函数
#
function WBICof(loglik)
    local E2nLn = makefunc_E2nLn(loglik)
    local n = size(loglik, 2)
    local WBIC = E2nLn(1/log(n))
    return WBIC
end

# 汎化損失を計算する函数
#
function GeneralizationLossof(pdftrue, pdfpred; xmin=-10.0, xmax=10.0)
    local f(x) = -pdftrue(x)*log(pdfpred(x))
    return quadgk(f, xmin, xmax)
end

# Shannon情報量を計算する函数
#
ShannonInformationof(pdftrue; xmin=-10.0, xmax=10.0) = GeneralizationLossof(pdftrue, pdftrue, xmin=xmin, xmax=xmax)[1]
ShannonInformationof(dist::Distribution, n) = entropy(dist)

# 最尤法を実行して AIC と BIC を計算する函数
#
# lpdf(w,y) = log p(y|w)
#
# w = link_model(z)   ←実ベクトル z 全体をパラメーター w の空間内に移す函数
# z = unlink_model(w) ←link_model(z)の逆函数
# これらは optimize() 函数の適用時にパラメーター w が定義域から外れることを防ぐための函数達
#
# (X[i], Y[i] はサンプル
#
# chain は loglik, chain = loglikchainof(sim) で作った chain
#
# optimize函数は1変数函数の場合には初期条件を与えて解を求める函数ではなくなるので、
# その場合には2変数函数に拡張してから使用している.
#
function AICandBICof(lpdf, link_model, unlink_model, Y, chain)
    local n = length(Y)
    local L = size(chain,1)
    local nparams = size(chain,2)
    local negloglik(z) = -@sum(lpdf(link_model(z), Y[i]), i, 1:n)
    local minnegloglik_chain, l
    minnegloglik_chain, l = findmin(negloglik(unlink_model(chain[l,:])) for l in 1:L)
    local o, minnegloglik, param_AIC
    if size(chain,2) == 1
        local f(z) = negloglik([z[1]]) + z[2]^2/eps()
        o = optimize(f, [unlink_model(chain[l,:])[1], 0.0])
        minnegloglik = o.minimum
        param_AIC = link_model([o.minimizer[1]])
    else
        o = optimize(negloglik, unlink_model(chain[l,:]))
        minnegloglik = o.minimum
        param_AIC = link_model(o.minimizer)
    end
    local T_AIC = 2.0*minnegloglik
    local V_AIC = 2.0*nparams
    local T_BIC = T_AIC
    local V_BIC = nparams*log(n)
    local AIC = T_AIC + V_AIC
    local BIC = T_BIC + V_BIC
    return AIC, T_AIC, V_AIC, 
        BIC, T_BIC, V_BIC, 
        param_AIC
end

function statsof(sim, Y; 
        dist_true=mixnormal(0.0, 0.0, 0.0), 
        dist_model=mixnormal, link_model=link_mixnormal, unlink_model=unlink_mixnormal)
    
    local n = length(Y)
    
    local lpdf(w, y) = logpdf(dist_model(w), y)
    
    local loglikchainof = makefunc_loglikchainof(lpdf, sort(keys(sim)), Y)
    local loglik, chain 
    loglik, chain = loglikchainof(sim)
    
    local WAIC, T_WAIC, V_WAIC
    WAIC, T_WAIC, V_WAIC = WAICof(loglik)
    local LOOCV = LOOCVof(loglik)
    
    local WBIC = WBICof(loglik)
    local FreeEnergy = FreeEnergyof(loglik)[1]
    
    local param_Bayes = vec(mean(chain, 1))
    local pred_Bayes = makefunc_pdfpred(lpdf, chain)
    
    local GeneralizationLoss = 2*n*GeneralizationLossof(x->pdf(dist_true,x), pred_Bayes)[1]
    
    local AIC, T_AIC, V_AIC, BIC, T_BIC, V_BIC, param_MLE
    AIC, T_AIC, V_AIC, BIC, T_BIC, V_BIC, param_MLE = AICandBICof(lpdf, link_model, unlink_model, Y, chain)
    
    local pred_MLE(y) = exp(lpdf(param_MLE, y))
    return WAIC, T_WAIC, V_WAIC, LOOCV, GeneralizationLoss,
        WBIC, FreeEnergy,
        param_Bayes, pred_Bayes, 
        AIC, T_AIC, V_AIC, 
        BIC, T_BIC, V_BIC,
        param_MLE, pred_MLE
end

function show_all_results(dist_true, Y, sim; statsfunc=statsof, 
        dist_model=mixnormal, link_model=link_mixnormal, unlink_model=link_mixnormal,
        figsize=(6,4.2), xmin=-4.0, xmax=6.0)
    WAIC, T_WAIC, V_WAIC, LOOCV, GeneralizationLoss,
    WBIC, FreeEnergy,
    param_Bayes, pred_Bayes, 
    AIC, T_AIC, V_AIC, 
    BIC, T_BIC, V_BIC,
    param_MLE, pred_MLE = statsfunc(sim, Y, dist_true=dist_true, 
        dist_model=dist_model, link_model=link_model, unlink_model=unlink_model)
    
    n = length(Y)
    println("\n=== Estimates by $dist_model  (n = $n) ===")
    @show param_Bayes
    @show param_MLE
    
    println("--- Information Criterions")
    println("* AIC     = $AIC = $T_AIC + $V_AIC")
    println("* GenLoss = $GeneralizationLoss")
    println("* WAIC    = $WAIC = $T_WAIC + $V_WAIC")
    println("* LOOCV   = $LOOCV")
    println("---")
    println("* BIC        = $BIC = $T_BIC + $V_BIC")
    println("* FreeEnergy = $FreeEnergy")
    println("* WBIC       = $WBIC")

    println("="^78 * "\n")
    
    sleep(0.1)
    plt.figure(figsize=figsize)
    kde_sample = makefunc_pdfkde(Y)
    x = linspace(xmin, xmax, 201)
    plt.plot(x, pdf.(dist_true, x), label="true distribution")
    plt.scatter(Y, kde_sample.(Y), label="sample", s=10, color="k", alpha=0.5)
    plt.plot(x, kde_sample.(x), label="KDE of sample",   color="k", alpha=0.5)
    plt.plot(x, pred_Bayes.(x), label="Baysian predictive", ls="--")
    plt.plot(x, pred_MLE.(x),   label="MLE predictive",     ls="-.")
    plt.xlabel("x")
    plt.ylabel("probability density")
    plt.grid(ls=":")
    plt.legend(fontsize=8)
    plt.title("Estimates by $dist_model: n = $(length(Y))")
end

function plotsample(dist_true, Y; figsize=(6,4.2), xmin=-4.0, xmax=6.0)
    sleep(0.1)
    plt.figure(figsize=figsize)
    kde_sample = makefunc_pdfkde(Y)
    x = linspace(xmin, xmax, 201)
    plt.plot(x, pdf.(dist_true, x), label="true dist.")
    plt.scatter(Y, kde_sample.(Y), label="sample", s=10, color="k", alpha=0.5)
    plt.plot(x, kde_sample.(x), label="KDE of sample",   color="k", alpha=0.5)
    plt.xlabel("x")
    plt.ylabel("probability density")
    plt.grid(ls=":")
    plt.legend(fontsize=8)
    plt.title("Sample size n = $(length(Y))")
end

# dist_true 分布に従う乱数で生成したサンプルの配列 Y を与えると
# 正規分布モデルの最尤法で推定して、AICなどを返す函数
#
function fit_Normal(dist_true, Y)
    local n = length(Y)
    local d = fit(Normal,Y)
    local mu = d.μ
    local sigma = d.σ
    local pred_Normal(y) = pdf(d,y)
    local T_Normal = -2*sum(logpdf.(d,Y))
    local V_Normal = 4.0
    local AIC_Normal = T_Normal + V_Normal
    local TB_Normal = T_Normal
    local VB_Normal = 2.0*log(n)
    local BIC_Normal = TB_Normal + VB_Normal
    local f(y) = -pdf(dist_true, y)*logpdf(d, y)
    local GL_Normal = 2*n*quadgk(f, -10, 10)[1]
    return mu, sigma, pred_Normal, 
        AIC_Normal, T_Normal, V_Normal, 
        BIC_Normal, TB_Normal, VB_Normal, 
        GL_Normal
end

# グラフをプロットする範囲が xmin から xmax まで
#
function show_fit_Normal(dist_true, Y; figsize=(6,4.2), xmin=-4.0, xmax=6.0)
    mu, sigma, pred_Normal, 
    AIC_Normal, T_Normal, V_Normal, 
    BIC_Normal, TB_Normal, VB_Normal, 
    GL_Normal = fit_Normal(dist_true, Y)
    println("--- Normal Fitting")
    println("* μ = $mu")
    println("* σ = $sigma")
    println("* GenLoss = $GL_Normal")
    println("* AIC     = $AIC_Normal = $T_Normal + $V_Normal")
    println("* BIC     = $BIC_Normal = $TB_Normal + $VB_Normal")
    
    sleep(0.1)
    plt.figure(figsize=figsize)
    kde_sample = makefunc_pdfkde(Y)
    x = linspace(xmin, xmax, 201)
    plt.plot(x, pdf.(dist_true, x), label="true distribution")
    plt.scatter(Y, kde_sample.(Y), label="sample", s=10, color="k", alpha=0.5)
    plt.plot(x, kde_sample.(x), label="KDE of sample",   color="k", alpha=0.5)
    plt.plot(x, pred_Normal.(x), label="Normal predictive", ls="--")
    plt.xlabel("x")
    plt.ylabel("probability density")
    plt.grid(ls=":")
    plt.legend(fontsize=8)
    plt.title("Sample size n = $(length(Y))")
end

function InformationCriterions(dist_true, Y, sim; dist_model=mixnormal)
    local n = length(Y)
    local lpdf(w, y) = logpdf(dist_model(w), y)

    local loglikchainof = makefunc_loglikchainof(lpdf, sort(keys(sim)), Y)
    local loglik, chain 
    loglik, chain = loglikchainof(sim)

    local WAIC, T_WAIC, V_WAIC
    WAIC, T_WAIC, V_WAIC = WAICof(loglik)
    local LOOCV = LOOCVof(loglik)
    
    local pred_Bayes = makefunc_pdfpred(lpdf, chain)
    local GL = 2*n*GeneralizationLossof(x->pdf(dist_true,x), pred_Bayes)[1]
    
    local WBIC = WBICof(loglik)
    local FreeEnergy = FreeEnergyof(loglik)[1]

    return chain, WAIC, T_WAIC, V_WAIC, LOOCV, GL, WBIC, FreeEnergy
 end

# サンプルを生成する真の確率分布の定義

theta = 0.4
dist_true = Gamma(1/theta^2, theta)
@show mu_true = mean(dist_true)
@show sigma_true = std(dist_true)
dist_normal_fitting = Normal(mu_true, sigma_true)
sleep(0.1)
x = linspace(-1,8,100)
plt.figure(figsize=(5,3.5))
plt.plot(x, pdf.(dist_true, x), label="true distribution")
plt.plot(x, pdf.(dist_normal_fitting, x), label="normal fitting")
plt.grid(ls=":")
plt.legend()

# シミュレーション

seed = 4649
dataname = "Gamma128Sample$seed"
N = 2^7

##########

iterations_mixnormal = 500
burnin_mixnormal = 100
thin_mixnormal = 1
chains_mixnormal = 4

iterations_normal1 = 500
burnin_normal1 = 100
thin_normal1 = 1
chains_normal1 = 4

iterations_normal = 500
burnin_normal = 100
thin_normal = 1
chains_normal = 4

Sample = rand(dist_true, N)
Shannon = Array{Float64}(N)
T_true = Array{Float64}(N)

Chain_mixnormal = Array{Array{Float64,2}}(N)
WAIC_mixnormal = Array{Float64}(N)
T_mixnormal = Array{Float64}(N)
V_mixnormal = Array{Float64}(N)
LOOCV_mixnormal = Array{Float64}(N)
GL_mixnormal = Array{Float64}(N)
WBIC_mixnormal = Array{Float64}(N)
FreeEnergy_mixnormal = Array{Float64}(N)

Chain_normal1 = Array{Array{Float64,2}}(N)
WAIC_normal1 = Array{Float64}(N)
T_normal1 = Array{Float64}(N)
V_normal1 = Array{Float64}(N)
LOOCV_normal1 = Array{Float64}(N)
GL_normal1 = Array{Float64}(N)
WBIC_normal1 = Array{Float64}(N)
FreeEnergy_normal1 = Array{Float64}(N)

Chain_normal = Array{Array{Float64,2}}(N)
WAIC_normal = Array{Float64}(N)
T_normal = Array{Float64}(N)
V_normal = Array{Float64}(N)
LOOCV_normal = Array{Float64}(N)
GL_normal = Array{Float64}(N)
WBIC_normal = Array{Float64}(N)
FreeEnergy_normal = Array{Float64}(N)

@time for i in 1:N
    n = i
    Y = Sample[1:n]
    Shannon[i] = 2*n*entropy(dist_true)
    T_true[i] = -2*@sum(logpdf(dist_true, Y[k]), k, 1:n)

    # mixnormal
    #
    model_mixnormal, data_mixnormal, inits_mixnormal = sample2model_mixnormal(Y, chains=chains_mixnormal)
    sim_mixnormal = mcmc(model_mixnormal, data_mixnormal, inits_mixnormal,
        iterations_mixnormal, burnin=burnin_mixnormal, thin=thin_mixnormal, chains=chains_mixnormal, verbose=false)

    Chain_mixnormal[i], WAIC_mixnormal[i], T_mixnormal[i], V_mixnormal[i], 
    LOOCV_mixnormal[i], GL_mixnormal[i], WBIC_mixnormal[i], FreeEnergy_mixnormal[i] = 
    InformationCriterions(dist_true, Y, sim_mixnormal, dist_model=mixnormal)
    
    # normal1
    #
    model_normal1, data_normal1, inits_normal1 = sample2model_normal1(Y, chains=chains_normal1)
    sim_normal1 = mcmc(model_normal1, data_normal1, inits_normal1,
        iterations_normal1, burnin=burnin_normal1, thin=thin_normal1, chains=chains_normal1, verbose=false)

    Chain_normal1[i], WAIC_normal1[i], T_normal1[i], V_normal1[i], 
    LOOCV_normal1[i], GL_normal1[i], WBIC_normal1[i], FreeEnergy_normal1[i] = 
    InformationCriterions(dist_true, Y, sim_normal1, dist_model=normal1)

    # normal
    #
    model_normal, data_normal, inits_normal = sample2model_normal(Y, chains=chains_normal)
    sim_normal = mcmc(model_normal, data_normal, inits_normal,
        iterations_normal, burnin=burnin_normal, thin=thin_normal, chains=chains_normal, verbose=false)

    Chain_normal[i], WAIC_normal[i], T_normal[i], V_normal[i], 
    LOOCV_normal[i], GL_normal[i], WBIC_normal[i], FreeEnergy_normal[i] = 
    InformationCriterions(dist_true, Y, sim_normal, dist_model=normal)
    
    print("$i ")
    i == N && println()
end

#########

varnames = [
    "seed", "dist_true", "N", "Shannon",
    "iterations_mixnormal", "burnin_mixnormal", "thin_mixnormal", "chains_mixnormal",
    "iterations_normal1",   "burnin_normal1",   "thin_normal1",   "chains_normal1",
    "iterations_normal",    "burnin_normal",    "thin_normal",    "chains_normal",
    "Sample", "T_true", 
    "Chain_mixnormal", "WAIC_mixnormal", "T_mixnormal",    "V_mixnormal", 
    "LOOCV_mixnormal", "GL_mixnormal",   "WBIC_mixnormal", "FreeEnergy_mixnormal", 
    "Chain_normal1",   "WAIC_normal1",   "T_normal1",      "V_normal1", 
    "LOOCV_normal1",   "GL_normal1",     "WBIC_normal1",   "FreeEnergy_normal1", 
    "Chain_normal",    "WAIC_normal",    "T_normal",       "V_normal", 
    "LOOCV_normal",    "GL_normal",      "WBIC_normal",    "FreeEnergy_normal", 
]
eval(parse("$dataname = Dict()"))
for v in sort(varnames)
    eval(parse("$dataname[\"$v\"] = $v"))
end

@show jld2file = "$dataname.jld2"
eval(parse("save(jld2file, $dataname)"))
data = load("$dataname.jld2")
@show keys(data);

# シミュレーション

seed = 12345
dataname = "Gamma128Sample$seed"
N = 2^7

##########

iterations_mixnormal = 500
burnin_mixnormal = 100
thin_mixnormal = 1
chains_mixnormal = 4

iterations_normal1 = 500
burnin_normal1 = 100
thin_normal1 = 1
chains_normal1 = 4

iterations_normal = 500
burnin_normal = 100
thin_normal = 1
chains_normal = 4

Sample = rand(dist_true, N)
Shannon = Array{Float64}(N)
T_true = Array{Float64}(N)

Chain_mixnormal = Array{Array{Float64,2}}(N)
WAIC_mixnormal = Array{Float64}(N)
T_mixnormal = Array{Float64}(N)
V_mixnormal = Array{Float64}(N)
LOOCV_mixnormal = Array{Float64}(N)
GL_mixnormal = Array{Float64}(N)
WBIC_mixnormal = Array{Float64}(N)
FreeEnergy_mixnormal = Array{Float64}(N)

Chain_normal1 = Array{Array{Float64,2}}(N)
WAIC_normal1 = Array{Float64}(N)
T_normal1 = Array{Float64}(N)
V_normal1 = Array{Float64}(N)
LOOCV_normal1 = Array{Float64}(N)
GL_normal1 = Array{Float64}(N)
WBIC_normal1 = Array{Float64}(N)
FreeEnergy_normal1 = Array{Float64}(N)

Chain_normal = Array{Array{Float64,2}}(N)
WAIC_normal = Array{Float64}(N)
T_normal = Array{Float64}(N)
V_normal = Array{Float64}(N)
LOOCV_normal = Array{Float64}(N)
GL_normal = Array{Float64}(N)
WBIC_normal = Array{Float64}(N)
FreeEnergy_normal = Array{Float64}(N)

@time for i in 1:N
    n = i
    Y = Sample[1:n]
    Shannon[i] = 2*n*entropy(dist_true)
    T_true[i] = -2*@sum(logpdf(dist_true, Y[k]), k, 1:n)

    # mixnormal
    #
    model_mixnormal, data_mixnormal, inits_mixnormal = sample2model_mixnormal(Y, chains=chains_mixnormal)
    sim_mixnormal = mcmc(model_mixnormal, data_mixnormal, inits_mixnormal,
        iterations_mixnormal, burnin=burnin_mixnormal, thin=thin_mixnormal, chains=chains_mixnormal, verbose=false)

    Chain_mixnormal[i], WAIC_mixnormal[i], T_mixnormal[i], V_mixnormal[i], 
    LOOCV_mixnormal[i], GL_mixnormal[i], WBIC_mixnormal[i], FreeEnergy_mixnormal[i] = 
    InformationCriterions(dist_true, Y, sim_mixnormal, dist_model=mixnormal)
    
    # normal1
    #
    model_normal1, data_normal1, inits_normal1 = sample2model_normal1(Y, chains=chains_normal1)
    sim_normal1 = mcmc(model_normal1, data_normal1, inits_normal1,
        iterations_normal1, burnin=burnin_normal1, thin=thin_normal1, chains=chains_normal1, verbose=false)

    Chain_normal1[i], WAIC_normal1[i], T_normal1[i], V_normal1[i], 
    LOOCV_normal1[i], GL_normal1[i], WBIC_normal1[i], FreeEnergy_normal1[i] = 
    InformationCriterions(dist_true, Y, sim_normal1, dist_model=normal1)

    # normal
    #
    model_normal, data_normal, inits_normal = sample2model_normal(Y, chains=chains_normal)
    sim_normal = mcmc(model_normal, data_normal, inits_normal,
        iterations_normal, burnin=burnin_normal, thin=thin_normal, chains=chains_normal, verbose=false)

    Chain_normal[i], WAIC_normal[i], T_normal[i], V_normal[i], 
    LOOCV_normal[i], GL_normal[i], WBIC_normal[i], FreeEnergy_normal[i] = 
    InformationCriterions(dist_true, Y, sim_normal, dist_model=normal)
    
    print("$i ")
    i == N && println()
end

#########

varnames = [
    "seed", "dist_true", "N", "Shannon",
    "iterations_mixnormal", "burnin_mixnormal", "thin_mixnormal", "chains_mixnormal",
    "iterations_normal1",   "burnin_normal1",   "thin_normal1",   "chains_normal1",
    "iterations_normal",    "burnin_normal",    "thin_normal",    "chains_normal",
    "Sample", "T_true", 
    "Chain_mixnormal", "WAIC_mixnormal", "T_mixnormal",    "V_mixnormal", 
    "LOOCV_mixnormal", "GL_mixnormal",   "WBIC_mixnormal", "FreeEnergy_mixnormal", 
    "Chain_normal1",   "WAIC_normal1",   "T_normal1",      "V_normal1", 
    "LOOCV_normal1",   "GL_normal1",     "WBIC_normal1",   "FreeEnergy_normal1", 
    "Chain_normal",    "WAIC_normal",    "T_normal",       "V_normal", 
    "LOOCV_normal",    "GL_normal",      "WBIC_normal",    "FreeEnergy_normal", 
]
eval(parse("$dataname = Dict()"))
for v in sort(varnames)
    eval(parse("$dataname[\"$v\"] = $v"))
end

@show jld2file = "$dataname.jld2"
eval(parse("save(jld2file, $dataname)"))
data = load("$dataname.jld2")
@show keys(data);