Table of Contents
混合正規分布モデルと正規分布モデルの各種情報量規準の比較¶
黒木玄
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
からダウンロードできる. 分析結果は以下の場所で閲覧できる:
n=8 http://nbviewer.jupyter.org/gist/genkuroki/cee2962916ec9269c194468a3e5cf288
n=32 http://nbviewer.jupyter.org/gist/genkuroki/41dd5005f143bda6cc678a243bca41c0
n=128 http://nbviewer.jupyter.org/gist/genkuroki/d2978887f61287d54e58baa700d2c002
警告・注意 (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
==============================================================================
