混合正規分布モデルと正規分布モデルの各種情報量規準の比較

黒木玄

2017-11-16~2018-01-13

正規分布の確率密度函数を次のように書くことにする:

$$\displaystyle N(\mu,\sigma) = \frac{\exp(-(x-\mu)^2/(2\sigma^2))}{\sqrt{2\pi\sigma^2}}. $$

標準正規分布 $N(0,1)$ で生成したサンプルについて, 分散を1に固定した山が2つの混合正規分布モデル mixnormal

$$\displaystyle y \sim (1-a)N(b,1) + a N(c,1) $$

と分散を1に固定した正規分布モデル normal1

$$\displaystyle y \sim N(\mu,1) $$

と分散も動かすパラメーター数2の正規分布モデル normal

$$\displaystyle y \sim N(\mu,\sigma) $$

による推定を比較する.

このノートで作成したデータファイルは

https://genkuroki.github.io/documents/Jupyter/#2017-11-16

からダウンロードできる. 分析結果は以下の場所で閲覧できる:

警告・注意 (2018-01-13) このノートブックにおけるWBICの数値計算法は適切ではない. なぜならば, このノートブックではWBICを逆温度β=1の場合に純粋数学的には成立しているWBICの公式を利用しているからである. 実際にやってみてわかったことだが(よくよく考えると当たり前のことだが), 逆温度β=1/log(n)における事後分布平均で定義されるWBICを逆温度β=1での事後分布のサンプルを使って近似計算すると誤差が大きくなる.

Julia言語のMamba.jlパッケージで逆温度βでのベイズ推定を行う方法については

http://nbviewer.jupyter.org/gist/genkuroki/8e97ef72ec01857564bd1a68e6e63d64

を参照せよ. 単純な逆温度βの正規分布モデルの場合にWBICのexact formulaを使った実験については

http://nbviewer.jupyter.org/gist/genkuroki/8a342d0b7b249e279dd8ad6ae283c5db

を参照せよ.

In [1]:
versioninfo()
Julia Version 0.6.0
Commit 903644385b* (2017-06-19 13:05 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

パッケージの読み込みと諸定義

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
 60.473134 seconds (15.80 M allocations: 1.103 GiB, 0.79% gc time)
Out[2]:
@sum (macro with 1 method)

MambaパッケージでMCMCを実行するための函数

In [3]:
@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
    )
    data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :a0 => a0,
        :b0 => b0,
        :c0 => c0,
    )
    
    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),
    )
    scheme = [
        NUTS([:a, :b, :c])
    ]
    setsamplers!(model, scheme)
    
    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
Out[3]:
sample2model_mixnormal (generic function with 1 method)
In [4]:
@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
    )
    data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
        :sigma0 => sigma0,
    )
    
    model = Model(
        y = Stochastic(1, (mu, sigma) -> dist_model(mu, sigma), false),
        mu = Stochastic(() -> prior_mu, true),
        sigma = Stochastic(() -> prior_sigma, true),
    )
    scheme = [
        NUTS([:mu, :sigma])
    ]
    setsamplers!(model, scheme)
    
    inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
            :sigma => data[:sigma0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end
Out[4]:
sample2model_normal (generic function with 1 method)
In [5]:
@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
    )
    data = Dict{Symbol, Any}(
        :Y => Y,
        :n => length(Y),
        :mu0 => mu0,
    )
    
    model = Model(
        y = Stochastic(1, mu -> dist_model(mu), false),
        mu = Stochastic(() -> prior_mu, true),
    )
    scheme = [
        NUTS([:mu])
    ]
    setsamplers!(model, scheme)
    
    inits = [
        Dict{Symbol, Any}(
            :y => data[:Y],
            :mu => data[:mu0],
        )
        for k in 1:chains
    ]
    return model, data, inits
end
Out[5]:
sample2model_normal1 (generic function with 1 method)

Mambaによるシミュレーション結果のまとめの表示

In [6]:
## 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)
        xmin = quantile(vec(sim[:,var,:].value), 0.005)
        xmax = quantile(vec(sim[:,var,:].value), 0.995)
        for k in sim.chains
            chain = sim[:,var,:].value[:,1,k]
            pdfkde = makefunc_pdfkde(chain)
            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
    m = 0
    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))

            u = keys_sim[i]
            v = keys_sim[j]
            X = vec(sim[:,u,:].value)
            Y = vec(sim[:,v,:].value)
            pdfkde = makefunc_pdfkde(X,Y)
            xmin = quantile(X, 0.005)
            xmax = quantile(X, 0.995)
            ymin = quantile(Y, 0.005)
            ymax = quantile(Y, 0.995)
            x = linspace(xmin, xmax, 200)
            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
Out[6]:
pyplot_contour (generic function with 1 method)

WAICなどを計算するための函数

以下において, 各種情報量規準の定義と計算の仕方を説明する.

情報量規準の表示のスケールには, Kullback-Leibler情報量のスケールと対数尤度比によるカイ二乗検定のスケールの2種類があるが, 以下では伝統的な後者のスケールに合わせて説明する. 対数尤度比によるカイ二乗検定のスケールからKL情報量のスケールに移るためには, サンプルサイズの2倍で割ればよい.

$X_1,\ldots,X_n$ はサイズ $n$ のサンプルであり, $d$ は確率モデル $p(x|w)$ のパラメーターの数であり($w=(w_1,\ldots,w_d$)), $\hat{w}$ は最大尤度を与えるパラメーターであるとする:

$$\displaystyle \min_w \left(-\sum_{i=1}^n \log p(X_k|w)\right) = -\sum_{i=1}^n \log p(X_k|\hat{w}). $$

最尤法の予測分布:

$$\displaystyle p^*_\mathrm{MLE}(x) = p(x|\hat{w}). $$

AICの計算の仕方:

$$\displaystyle \mathrm{AIC} := 2\left(-\sum_{k=1}^n\log p(X_k|\hat{w})\right) + 2d. $$

BICの計算の仕方:

$$\displaystyle \mathrm{BIC} := 2\left(-\sum_{k=1}^n\log p(X_k|\hat{w})\right) + d\log n. $$

MCMCの結果得られる事後分布のサンプルを $w_1,w_2,\ldots,w_L$ と書く.

事後予測分布とWAICとLOOCV(一個抜き交差検証)とWBICは対数尤度の配列

$$ \{\,\log p(X_k|w_l)\mid k=1,2,\ldots,n,\; l=1,2,\ldots,L\,\} $$

のみから近似計算可能である. 事後予測分布, WAIC, LOOCV, WBIC を計算するための最初の作業はこの配列をMCMCの結果から抽出することである. その配列を抽出できてさえしまえば, それらは以下の公式によって機械的に計算される.

Bayes推定の事後予測分布:

$$\displaystyle p^*_\mathrm{Bayesian}(x) := (\mathrm{posterior\ mean\ of}\ p(x|w)) \approx \left(\mathrm{mean\ of}\ \{p(x|w_l)\}_{l=1}^L\right) = \left(\mathrm{mean\ of}\ \{\exp(\log p(x|w_l))\}_{l=1}^L\right). $$

WAICの計算の仕方:

$$\displaystyle \mathrm{WAIC} := T_{\mathrm{WAIC}} + V_{\mathrm{WAIC}}. $$

ここで $$\begin{aligned} & T_{\mathrm{WAIC}} := -2\sum_{k=1}^n \log p^*(X_k) \approx -2\sum_{k=1}^n\log\left(\mathrm{mean\ of}\ \{\exp(\log p(X_k|w_l))\}_{l=1}^L\right), \\ & V_{\mathrm{WAIC}} := 2\sum_{k=1}^n \left(\mathrm{posterior\ variance\ of}\ \log p(X_k|w)\right) \approx 2\sum_{k=1}^n \left(\mathrm{variance\ of}\ \{\log p(X_k|w_l)\}_{l=1}^L\right). \end{aligned}$$

LOOCV(一個抜き交差検証)の素朴な計算の仕方

$$\begin{aligned} \mathrm{LOOCV} &:= -2\sum_{k=1}^n \log E^{(k)}_w[p(X_k|w)] = 2\sum_{k=1}^n \log E_w\left[\frac{1}{p(X_k|w)}\right] \\ & \approx 2\sum_{k=1}^n \log\left(\mathrm{mean\ of}\ \{\exp(-\log p(X_k|w_l))\}_{l=1}^L\right). \end{aligned}$$

ここで

$$\begin{aligned} & Z_n := \int \prod_{j=1}^n p(X_j|w)\;\varphi(w)\,dw, \\ & E_w[f(w)] := \frac{1}{Z_n}\int f(w)\prod_{j=1}^n p(X_j|w)\;\varphi(w)\,dw \approx \left(\mathrm{mean\ of}\ \{f(w_l)\}_{l=1}^L\right), \\ & Z^{(k)}_n := \int \frac{\prod_{j=1}^n p(X_j|w)}{p(X_k|w)}\varphi(w)\,dw = Z_n E_w\left[\frac{1}{p(X_k|w)}\right], \\ & E^{(k)}_w[f(w)] := \frac{1}{Z^{(k)}_n}\int f(w)\frac{\prod_{j=1}^n p(X_j|w)}{p(X_k|w)}\varphi(w)\,dw = \frac{Z_n}{Z^{(k)}_n}E_w\left[\frac{f(w)}{p(X_k|w)}\right] = \frac{E_w[f(w)/p(X_k|w)]}{E_w[1/p(X_k|w)]}, \\ & \log E^{(k)}_w[p(X_k|w)] = \log\frac{E_w[1]}{E_w[1/p(X_k|w)]} = -\log E_w\left[\frac{1}{p(X_k|w)}\right] \\ & \qquad \approx -\log\left(\mathrm{mean\ of}\ \left\{\frac{1}{p(X_k|w_l)}\right\}_{l=1}^L\right) = -\log\left(\mathrm{mean\ of}\ \{\exp(-\log p(X_k|w_l))\}_{l=1}^L\right). \end{aligned}$$

WBICの計算の仕方: $\beta=1/\log n$ とおき,

$$\displaystyle \mathrm{WBIC} := 2E^{\beta}_w\left[nL_n(w)\right] = 2\frac{N_{\mathrm{WBIC}}}{D_{\mathrm{WBIC}}}. $$

ここで,

$$\begin{aligned} nL_n(w) &:= -\sum_{k=1}^n \log p(X_k|w), \quad nL_n(w_l) = -\sum_{k=1}^n \log p(X_k|w_l), \\ Z_n(\beta) &:= \int \left(\prod_{k=1}^n p(X_k|w)\right)^\beta\;\varphi(w)\,dw, \quad Z_n(\beta)E^{\beta}_w\left[f(w)\right] := \int f(w)\left(\prod_{k=1}^n p(X_k|w)\right)^\beta\;\varphi(w)\,dw \\ N_{\mathrm{WBIC}} &:= \frac{Z_n(\beta)E^{\beta}_w\left[nL_n(w)\right]}{Z_n(1)} = \frac{1}{Z_n(1)} \int nL_n(w)\left(\prod_{k=1}^n p(X_k|w)\right)^{\beta-1} \prod_{k=1}^n p(X_k|w)\;\varphi(w)\,dw \\ & \approx \left(\mathrm{mean\ of}\ \left\{nL_n(w_l) \exp(-(\beta-1)nL_n(w_l))\right\}_{l=1}^L\right), \\ D_{\mathrm{WBIC}} &:= \frac{Z_n(\beta)}{Z_n(1)} = \frac{1}{Z_n(1)}\int \left(\prod_{k=1}^n p(X_k|w)\right)^{\beta-1} \prod_{k=1}^n p(X_k|w)\;\varphi(w)\,dw \\ & \approx \left(\mathrm{mean\ of}\ \left\{\exp(-(\beta-1)nL_n(w_l))\right\}_{l=1}^L\right). \end{aligned}$$

WBICは逆温度 $1$ の事後分布のサンプル上での対数尤度函数の $-1$ 倍の値達 $nL_n(w_l)$ のみから近似計算可能である.

WBICの近似計算のためには逆温度 $\beta$ の事後分布のサンプルを構成する必要がないことに注意せよ!

文献: 以上のWBICの計算法では,

http://www.jmlr.org/papers/volume14/watanabe13a/watanabe13a.pdf p.877 (p.11 in file)

の式(20)を $\beta_1=1$, $\beta_2=\beta$ の場合に適用した結果を使った.

2017-11-10 メモ: このノートブックにおけるWBICの計算はMCMCの結果に大きく依存し、分散が大きい。もしかしたら、何か間違ったことをしているかもしれない.

2017-11-12 追記 (自由エネルギーのより正確な計算法): WBIC導出の第一段階は, 自由エネルギー $F_n(\beta)$ の逆温度 $\beta$ による微分が $nL_n(w)$ の逆温度 $\beta$ の事後分布に関する平均になり, $F_n(0)=0$ なので

$$\displaystyle F_n(1) = \int_0^1 E^\beta_w[nL_n(w)] \,d\beta \tag{$*$} $$

となることを示すことである(これはとても易しい). WBIC導出の第二段階は右辺の積分が $\beta=1/\log n$ のときの $E^\beta[nL_n(w)]$ で近似できることを示すことである(こちらは少々非自明である).

($*$)の導出: $$ \begin{aligned} nL_n(w) &:= -\sum_{k=1}^n \log p(X_k|w), \\ Z_n(\beta) &:= \int \left(\prod_{k=1}^n p(X_k|w)\right)^\beta \varphi(w)\,dw = \int\exp\left(-\beta nL_n(w)\right)\varphi(w)\,dw, \\ Z_n'(\beta) &= -\int nL_n(w)\exp\left(-\beta nL_n(w)\right)\varphi(w)\,dw, \\ E^\beta_w[f(w)] &:= \frac{1}{Z_n(\beta)}\int f(w)\left(\prod_{k=1}^n p(X_k|w)\right)^\beta \varphi(w)\,dw = \frac{\int\exp\left(-\beta nL_n(w)\right)\varphi(w)\,dw}{Z_n(\beta)}, \\ F_n(\beta) &:= -\log Z_n(\beta), \\ F_n(0) &= -\log Z_n(0) = -\log\int\varphi(w)\,dw = -\log 1 = 0, \\ F_n'(\beta) &= \frac{-Z_n'(\beta)}{Z_n(\beta)} = \frac{\int nL_n(w)\exp\left(-\beta nL_n(w)\right)\varphi(w)\,dw}{Z_n(\beta)} = E^\beta_w[nL_n(w)]. \\ \therefore \quad F_n(1) &= \int_0^1 E^\beta_w[nL_n(w)]\,d\beta. \end{aligned} $$

$E^\beta_w[nL_n(w)]$ を $\beta$ で数値積分すれば自由エネルギー $F_n(1)$ を数値計算できる. 数値積分の分だけ必要な計算量は増えるが, 被積分函数は $\beta$ の1変数函数なので, モデルのパラメーターに関する多変数函数の積分をするよりも必要な計算量は圧倒的に少なくてすむ.

$E^\beta_w[nL_n(w)]$ は次のように表せる:

$$\displaystyle E^{\beta}_w\left[nL_n(w)\right] = \frac{N(\beta)}{D(\beta)}. $$

ここで,

$$\begin{aligned} N(\beta) &= \frac{Z_n(\beta)E^{\beta}_w\left[nL_n(w)\right]}{Z_n(1)} = \frac{1}{Z_n(1)} \int nL_n(w)\left(\prod_{k=1}^n p(X_k|w)\right)^{\beta-1} \prod_{k=1}^n p(X_k|w)\;\varphi(w)\,dw \\ & \approx \left(\mathrm{mean\ of}\ \left\{nL_n(w_l) \exp(-(\beta-1)nL_n(w_l))\right\}_{l=1}^L\right), \\ D(\beta) &= \frac{Z_n(\beta)}{Z_n(1)} = \frac{1}{Z_n(1)} \int \left(\prod_{k=1}^n p(X_k|w)\right)^{\beta-1} \prod_{k=1}^n p(X_k|w)\;\varphi(w)\,dw \\ & \approx \left(\mathrm{mean\ of}\ \left\{\exp(-(\beta-1)nL_n(w_l))\right\}_{l=1}^L\right). \end{aligned}$$

注意: 伝統的なAICのスケールに合わせるためには

$$\displaystyle 2F_n(1) = 2\int_0^1 E^{\beta}_w\left[nL_n(w)\right]\,d\beta = 2\int_0^1\frac{N(\beta)}{D(\beta)}\,d\beta $$

を計算すればよい. 以下のセルでは以上の方法でこのスケールで自由エネルギーを計算する函数を定義している.

In [7]:
# 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,:] が抽出される
    #
    function loglikchainof(sim)
        val = sim[:, symbols, :].value
        chain = vcat((val[:,:,k] for k in 1:size(val,3))...)
        L = size(chain,1)
        n = length(Y)
        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)
    L = size(chain,1)
    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)
    L, n = size(loglik)
    T_n = -mean(log(mean(exp(loglik[l,i]) for l in 1:L)) for i in 1:n)
    V_n  = sum(var(loglik[:,i], corrected=false) for i in 1:n)
    WAIC = 2*n*T_n + 2*V_n
    WAIC, 2*n*T_n, 2*V_n
end

# loglik[l,i] からLOOCVを素朴に計算する函数
#
function LOOCVof(loglik)
    L, n = size(loglik)
    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)
    E2nLn = makefunc_E2nLn(loglik)
    F = quadgk(E2nLn, 0.0, 1.0)[1]
    return F, 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

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

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

# 最尤法を実行して 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)
    n = length(Y)
    L = size(chain,1)
    nparams = size(chain,2)
    negloglik(z) = -@sum(lpdf(link_model(z), Y[i]), i, 1:n)
    minnegloglik_chain, l
    minnegloglik_chain, l = findmin(negloglik(unlink_model(chain[l,:])) for l in 1:L)
    if size(chain,2) == 1
        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
    T_AIC = 2.0*minnegloglik
    V_AIC = 2.0*nparams
    T_BIC = T_AIC
    V_BIC = nparams*log(n)
    AIC = T_AIC + V_AIC
    BIC = T_BIC + V_BIC
    return AIC, T_AIC, V_AIC, 
        BIC, T_BIC, V_BIC, 
        param_AIC
end
Out[7]:
AICandBICof (generic function with 1 method)

情報をまとめて表示するための函数

In [8]:
function statsof(sim, Y; 
        dist_true=mixnormal(0.0, 0.0, 0.0), 
        dist_model=mixnormal, link_model=link_mixnormal, unlink_model=unlink_mixnormal)
    
    n = length(Y)
    
    lpdf(w, y) = logpdf(dist_model(w), y)
    
    loglikchainof = makefunc_loglikchainof(lpdf, sort(keys(sim)), Y)
    loglik, chain = loglikchainof(sim)
    
    WAIC, T_WAIC, V_WAIC = WAICof(loglik)
    LOOCV = LOOCVof(loglik)
    
    WBIC = WBICof(loglik)
    FreeEnergy = FreeEnergyof(loglik)[1]
    
    param_Bayes = vec(mean(chain, 1))
    pred_Bayes = makefunc_pdfpred(lpdf, chain)
    
    GeneralizationLoss = 2*n*GeneralizationLossof(x->pdf(dist_true,x), pred_Bayes)[1]
    
    AIC, T_AIC, V_AIC, BIC, T_BIC, V_BIC, param_MLE = AICandBICof(lpdf, link_model, unlink_model, Y, chain)
    
    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
Out[8]:
plotsample (generic function with 1 method)

単純な正規分布モデルを最尤法で解くための函数

次のセルの定義はこのノートでは使用しない.

In [9]:
# dist_true 分布に従う乱数で生成したサンプルの配列 Y を与えると
# 正規分布モデルの最尤法で推定して、AICなどを返す函数
#
function fit_Normal(dist_true, Y)
    n = length(Y)
    d = fit(Normal,Y)
    mu = d.μ
    sigma = d.σ
    pred_Normal(y) = pdf(d,y)
    T_Normal = -2*sum(logpdf.(d,Y))
    V_Normal = 4.0
    AIC_Normal = T_Normal + V_Normal
    TB_Normal = T_Normal
    VB_Normal = 2.0*log(n)
    BIC_Normal = TB_Normal + VB_Normal
    f(y) = -pdf(dist_true, y)*logpdf(d, y)
    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
Out[9]:
show_fit_Normal (generic function with 1 method)

一回だけの比較

In [10]:
show_sample = false

macro simonce(seed, n, dist_true) quote begin

srand($(esc(seed)))

dist_true = $(esc(dist_true))
n = $(esc(n))
Y = rand(dist_true, n)

println("dist_true = $dist_true")
println("n = $n")
show_sample && println("Y = $Y")
println()

iterations = 500
burnin = 100
thin = 1
chains = 4
model, data, inits = sample2model_mixnormal(Y, chains=chains)
@time sim_mix = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim_mix, show_describe=false, show_gelman=false)
show_all_results(dist_true, Y, sim_mix, xmax=4.5,
    dist_model=mixnormal, link_model=link_mixnormal, unlink_model=unlink_mixnormal)

iterations = 500
burnin = 100
thin = 1
chains = 4
model, data, inits = sample2model_normal1(Y, chains=chains)
@time sim_normal1 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim_normal1, show_describe=false, show_gelman=false)
show_all_results(dist_true, Y, sim_normal1, xmax=4.5, 
    dist_model=normal1, link_model=link_normal1, unlink_model=unlink_normal1)

iterations = 500
burnin = 100
thin = 1
chains = 4
model, data, inits = sample2model_normal(Y, chains=chains)
@time sim_normal = mcmc(model, data, inits, iterations, burnin=burnin, thin=thin, chains=chains, verbose=false)
showsummary(sim_normal, show_describe=false, show_gelman=false)
show_all_results(dist_true, Y, sim_normal, xmax=4.5, 
    dist_model=normal, link_model=link_normal, unlink_model=unlink_normal)

end end end; # macro
In [11]:
@simonce(4649, 2^5, Normal())
dist_true = Distributions.Normal{Float64}(μ=0.0, σ=1.0)
n = 32

 15.065452 seconds (2.95 M allocations: 163.227 MiB, 0.55% gc time)

=== Estimates by mixnormal  (n = 32) ===
param_Bayes = [0.486959, 0.233942, 0.261794]
param_MLE = [0.936807, -1.19038, 0.501329]
--- Information Criterions
* AIC     = 98.63118881738633 = 92.63118881738633 + 6.0
* GenLoss = 95.00645794903944
* WAIC    = 95.23777068273637 = 93.032919310703 + 2.204851372033376
* LOOCV   = 95.24412271029585
---
* BIC        = 103.02839652578551 = 92.63118881738633 + 10.39720770839918
* FreeEnergy = 95.53777055940968
* WBIC       = 96.1060899100275
==============================================================================

  6.376224 seconds (15.31 k allocations: 1.008 MiB)

=== Estimates by normal1  (n = 32) ===
param_Bayes = [0.378015]
param_MLE = [0.394422]
--- Information Criterions
* AIC     = 95.13734050276199 = 93.13734050276199 + 2.0
* GenLoss = 95.26561685218313
* WAIC    = 95.17694627294915 = 93.09311908752528 + 2.0838271854238655
* LOOCV   = 95.17759864932783
---
* BIC        = 96.60307640556171 = 93.13734050276199 + 3.4657359027997265
* FreeEnergy = 95.04674485334121
* WBIC       = 95.46863326332175
==============================================================================

  0.429267 seconds (22.90 k allocations: 1.495 MiB)

=== Estimates by normal  (n = 32) ===
param_Bayes = [0.387085, 1.09131]
param_MLE = [0.156444, 1.16935]
--- Information Criterions
* AIC     = 99.25302370284928 = 95.25302370284928 + 4.0
* GenLoss = 95.41341947518924
* WAIC    = 97.28749496835465 = 93.36300964008126 + 3.924485328273391
* LOOCV   = 97.34306189082106
---
* BIC        = 102.18449550844873 = 95.25302370284928 + 6.931471805599453
* FreeEnergy = 97.79522662090972
* WBIC       = 99.0407046290298
==============================================================================