要rstan,loo

はじめに

本ノートは，ベルヌイ確率モデルとベータ事前分布を仮定した単純例をもとに，AIC，WAIC，周辺尤度とその近似としてのBIC，WBICの意味を実感することを目的とする．

q<-0.7
n<-30
a<-1
b<-1

最初に例示データとして，確率\(q=0.7\)のベルヌイ分布からの\(n=30\)個のサンプルを想定する．

# ex<-rbinom(n,1,q)
ex<-c(1,1,1,1,1,1,0,1,1,1,1,1,0,1,0,1,0,0,0,0,1,1,0,1,0,0,0,1,1,1)

\[X_i \sim Bern(q)\]

ベイズモデルの事前分布としてはベータ分布を想定する． \[q\sim Beta(\alpha,\beta)\] とくに具体例においては\(\alpha=1\), \(\beta=1\)を仮定する．

plot(seq(0,1,0.1),dbeta(seq(0,1,0.1),a,b),
     xlab="q",ylab="prior dens.",type="l")

(1) AIC

一般的には，以下のことを想定する．

最尤推定量

\(X^n\)の実現値として\(x^n\)が与えられたとき，最尤推定量\(\hat{\theta}_{ML}\)は，確率モデルの\(\theta\)についての最大化問題として \[\hat{\theta}_{ML}=\mathop{\rm arg~max}_\theta p(x^n|\theta)\] と解くことができる．

予測分布

最尤推定量をもとにした真の分布についての予測分布は \[p^*(X)=p(X|\hat{\theta}_{ML})\] である．つまり，最尤推定量を確率モデルに入れ込んだものが，真の分布ではないかと推測しているのである．

ベルヌイ確率モデルを想定した具体例では，\(\bar{x}=\sum_i x_i/n\)のとき， \[p^*(X)=p(X|\bar{x})=(\bar{x})^X(1-\bar{x})^{1-X}\] である．

ML_pred_dist<-function(x,data) {
  (sum(data)/length(data))^x*(1-sum(data)/length(data))^(1-x)
}
ML_pred_dist(c(1,0),ex)
## [1] 0.6333333 0.3666667

汎化損失

一般に，予測分布\(p^*(X)\)の対数をとってマイナスしたもの（つまり予測分布における\(X\)の情報量）について，真の分布\(q(X)\)について平均をとったものを「汎化損失(Generalization loss)」という．つまり，汎化損失\(G_n\) は， \[G_n=-E_X[\log p^*(X)]\]

汎化損失\(G_n\)は，真の分布のエントロピーと真の分布と予測分布のカルバック=ライブラ(KL)距離に分解できる（渡辺 2012: 9-10）． \[\begin{align*} G_n&=-E_X[\log p^*(X)] =-\int \log p^*(x) q(x)dx \\ &=-\int \log q(x) q(x)dx+\int \log \frac{q(x)}{p^*(x)} q(x)dx \end{align*}\]

真の分布のエントロピーは変わらないとすると，真の分布と予測分布の平均的な遠さを表すKL距離が小さくなるほど，汎化損失は小さくなる．そこで，汎化損失を真の分布の予測の指標とすることができるだろう．

平均対数尤度

最尤推定のとき，とくに汎化損失は「平均対数尤度(Expected log likelihood)」の符号を逆転したものに一致する． \[G_n=-E_X[\log p^*(X)]=-E_X[\log p(X|\hat{\theta}_{ML})]=-\int \log p(x|\hat{\theta}_{ML}) q(x)dx\]

汎化損失が小さくなるほど，平均対数尤度は大きくなるので，平均対数尤度の大きさは予測の良さを示していると考えられるだろう．

具体例における平均対数尤度は，真の分布がわかっているので，以下のようになる． \[E_X[\log p^*(X)]=E_X[\log p(X|\bar{x})]=q \log \bar{x} +(1-q)\log(1-\bar{x}) \]

E_log_likelihood<-function(data) {
  q*log(sum(data)/length(data))+(1-q)*log(1-sum(data)/length(data))
}
E_log_likelihood(ex)
## [1] -0.6207215
## G_n
-E_log_likelihood(ex)
## [1] 0.6207215

この具体例では，真の分布がわかっているので汎化損失が一意的に求められたが，通常は未知である．この汎化損失をデータから推定することが目的となる．

経験損失

一般的に，予測分布\(p^*(x)\)の対数をとってマイナスしたものに，すでに予測分布の導出に用いたデータを入れて，平均をとったものを「経験損失(Training loss)」という． \[T_n=-\frac{1}{n}\sum_i \log p^*(x_i)\]

最大対数尤度

最尤推定において，経験損失は「最大対数尤度(maximum log likelihood)」に\(-1/n\)をかけて求められる．

\[T_n=-\frac{1}{n}\sum_i \log p(x_i|\hat{\theta}_{ML})\]

具体例における最大対数尤度は，以下のようになる． \[\sum_i \log p(x_i|\hat{\theta}_{ML})=\sum_i \left[x_i \log \bar{x} +(n-x_i)\log(1-\bar{x})\right]\]

maximum_log_likelihood<-function(data) {
  sum(log(ML_pred_dist(data,data)))
}
maximum_log_likelihood(ex)
## [1] -19.71473
## T_n
-(1/length(ex))*maximum_log_likelihood(ex)
## [1] 0.6571578

経験損失は，データが与えられたら計算できるので，これによって，汎化損失を推定することを考える．

AIC

予測分布導出に使ったデータと同じデータを再び入れる経験損失は，汎化損失から何らかの「偏り(bias)」が生じる．

つまり， \[b(x^n)=G_n-T_n=-E_X[\log p(X|\hat{\theta}_{ML})]+\frac{1}{n}\sum_i \log p(x_i|\hat{\theta}_{ML})\]

-E_log_likelihood(ex)+(1/length(ex))*maximum_log_likelihood(ex)
## [1] -0.03643625

このバイアスの期待値を求めてみよう．ここでは，\(\hat{\theta}_{ML}\)も\(X^n\sim q(X)\)によって確率的に動くことを強調するために，\(\hat{\theta}_{ML}(X^n)\)と表す．

\[E_{X^n}[b(X^n)]=E_{X^n}\left[-E_X[\log p(X|\hat{\theta}_{ML}(X^n))]+\frac{1}{n}\sum_i \log p(X_i|\hat{\theta}_{ML}(X^n))\right]\] ここではシミュレーションで求める．

bias<-function(data) {
  -E_log_likelihood(data)+(1/length(data))*maximum_log_likelihood(data)
}
bias_dist<-replicate(10000,bias(rbinom(n,1,q)))
hist(bias_dist,col="skyblue")

mean(bias_dist[!is.infinite(bias_dist)])
## [1] 0.03405859

だいたい0.034となった．これは，モデルのパラメータ数1をサンプルサイズ30でわったものと平均的に一致している．

じっさい，最大対数尤度と\(n\times\)平均対数尤度のバイアスは，「想定される確率モデルのなかに真の分布が含まれる」という仮定の下で漸近的に自由パラメータ数\(d\)に一致し（小西・北川 2004: 50-55），汎化損失と経験損失のバイアスは\(d/n\)に一致する．

経験損失にバイアス\(d/n\)を加えたものが「AIC(Akaike Information Criterion)」である（一般的な定義としてはさらに\(2n\)倍する）．

\[\mathrm{AIC}=-\frac{1}{n}\sum_i \log p(x_i|\hat{\theta}_{ML})+\frac{d}{n}\]

AIC<-function(data) {
  -(1/length(data))*maximum_log_likelihood(data)+1/length(data)
}
AIC(ex)
## [1] 0.6904911

確認のため，Rのglm関数でAICを計算する（glmのAICは\(2n\)倍したもの）．

g<-glm(ex ~ 1, family = binomial)
g$aic/(2*n)
## [1] 0.6904911

理論的に，次のことが成り立つ（渡辺 2012: 80）． \[E_{X^n}\left(\mathrm{AIC}\right)=E_{X^n}(G_n)+o(1/n)\]

では，\(AIC\)が平均的に汎化損失（平均対数尤度の符号逆転）のよい推定量になっているかを確認しよう．

## AIC
AIC_dist<-replicate(10000,AIC(rbinom(n,1,q)))

## -(expected log likelihood)
mll_dist<-replicate(10000,-E_log_likelihood(rbinom(n,1,q)))

## compare
boxplot(AIC_dist,mll_dist,
        names=c("AIC","-mll"),col=c("#FF00007F","#0000FF7F"))

c(mean(AIC_dist),mean(mll_dist[!is.infinite(mll_dist)]))
## [1] 0.6267180 0.6294989

一致かくにん！　よかった．

(2) WAIC

ベイズ推定では一般的に以下を想定する．

事後分布

\[p(\theta|X^n)=\frac{p(X^n|\theta)\varphi(\theta)}{\int p(X^n|\theta)\varphi(\theta)d\theta}=\frac{\prod_i p(X_i|\theta)\varphi(\theta)}{\int \prod_i p(X_i|\theta)\varphi(\theta)d\theta}\] \(X^n\)の実現値として\(x^n\)が与えられたときのベルヌイ確率モデルの事後分布は以下の通り．ただし，\(B(a,b)\)はベータ関数である． \[p(q|x^n)=\frac{\prod_i p(x_i|\theta)\varphi(\theta)}{\int \prod_i p(X_i|\theta)\varphi(\theta)d\theta}=\frac{q^{\alpha+\sum_i x_i-1}(1-q)^{\beta+n-\sum_i x_i-1}}{B(\alpha+\sum_i x_i,\beta+n-\sum_i x_i)} \] つまり，事後分布はベータ分布\(Beta(\alpha+\sum_i x_i,\beta+n-\sum_i x_i)\)である． \(\alpha=1\), \(\beta=1\)のとき，とくに事後分布は，\(Beta(1+\sum_i x_i,1+n-\sum_i x_i)\)．

plot(seq(0,1,0.01),dbeta(seq(0,1,0.01),a+sum(ex),b+length(ex)-sum(ex)),
     xlab="q",ylab="posterior dens.",type="l")

予測分布

ベイズ推定をもとにした真の分布についての予測分布は \[p^*(X)=E_\theta[p(X|\theta)]=\int p(X|\theta)p(\theta|X^n)d\theta\] である．つまり，確率モデルを事後分布によって期待値をとったものが，真の分布ではないかと推測しているのである．

ベルヌイ確率モデルを想定した具体例では， \[p^*(X)=\int_0^1 p(X|q)p(q|x^n)dq=\left(\frac{\alpha+\sum_i x_i}{\alpha+\beta+n}\right)^X\left(\frac{\beta+n-\sum_i x_i}{\alpha+\beta+n}\right)^{1-X}\] である．

Bayes_pred_dist<-function(x,data) {
  ((a+sum(data))/(a+b+length(data)))^x*
    ((b+length(data)-sum(data))/(a+b+length(data)))^(1-x)
}
Bayes_pred_dist(c(1,0),ex)
## [1] 0.625 0.375

汎化損失

ベイズ予測分布についての汎化損失\(G\)は， \[G_n=-E_X[\log p^*(X)]=-E_X[\log E_\theta[p(X|\theta)]]=-\int q(x)\log \int p(x|\theta)p(\theta|X^n)d\theta dx\] ベルヌイ確率モデルを想定した具体例では， \[G_n=-q\log\left(\frac{\alpha+\sum_i x_i}{\alpha+\beta+n}\right)-(1-q)\log\left(\frac{\beta+n-\sum_i x_i}{\alpha+\beta+n}\right)\]

gen_loss<-function(data) {
  -q*log((a+sum(data))/(a+b+length(data)))-
    (1-q)*log((b+length(data)-sum(data))/(a+b+length(data)))
}
gen_loss(ex)
## [1] 0.6232513

この具体例では，真の分布がわかっているので，あるデータに対する汎化損失が一意的に求められたが，通常は未知である．この汎化損失をデータから推定することが目的となる．

経験損失

ベイズ予測分布についての経験損失\(T_n\)は， \[T_n=-\frac{1}{n}\sum_i \log p^*(x_i)=-\frac{1}{n}\sum_i\log E_\theta[p(x_i|\theta)]\] ベルヌイ確率モデルを想定した具体例では， \[T_n=-\left(\frac{\sum_i x_i}{n}\right)\log\left(\frac{\alpha+\sum_i x_i}{\alpha+\beta+n}\right)-\left(1-\frac{\sum_i x_i}{n}\right)\log\left(\frac{\beta+n-\sum_i x_i}{\alpha+\beta+n}\right)\]

T_loss<-function(data) {
  -mean(log(Bayes_pred_dist(data,data)))
}
T_loss(ex)
## [1] 0.6573064

経験損失は，データが与えられたら計算できるので，経験損失で汎化損失を推定することを考えよう．

WAIC

予測分布導出に使ったデータと同じデータを再び入れる経験損失は，汎化損失から何らかの偏りが生じる．ここでは，汎化損失から経験損失がどれくらい離れているかを見よう．

つまり， \[\beta(x^n)=G_n(x^n)-T_n(x^n)\]

gen_loss(ex)-T_loss(ex)
## [1] -0.03405504

このバイアスの期待値を求めてみよう． \[E_{X^n}[\beta(X^n)]=E_{X^n}\left[G_n(X^n)-T_n(X^n)\right]\]

Bayes_bias<-function(data) {
  gen_loss(data)-T_loss(data)
}
Bayes_bias_dist<-replicate(10000,Bayes_bias(rbinom(n,1,q)))
hist(Bayes_bias_dist,col="skyblue")

mean(Bayes_bias_dist)
## [1] 0.03101397

だいたい0.0031となった．

汎化損失と経験損失の平均的なバイアスは，以下で定義される「汎関数分散(general function variance)」を\(n\)でわった\(V_n/n\)と漸近的に一致することが知られている（渡辺 2012: 118）．

\[V_n=\sum_i\left\{E_\theta[(\log p(x_i|\theta))^2]-E_\theta[\log p(x_i|\theta)]^2\right\}\] ベルヌイ確率モデルを想定した具体例では，以下のようになる（Mathematicaさんにやってもらった）． \[\begin{align*} V_n&=\sum_i\left\{\int (x_i \log q+(1-x_i) \log(1-q))^2p(q|x^n) dq -\left(\int (x_i \log q+(1-x_i) \log(1-q))p(q|x^n) dq\right)^2\right\} \\ &=\sum_i\left\{x_i \left[\psi^1(\alpha+\sum_j x_j)-\psi^1(\alpha+\beta+n)\right]+(1-x_i) \left[\psi^1(\beta+n-\sum_j x_j)-\psi^1(\alpha+\beta+n)\right]\right\} \\ &=\left(\sum_i\ x_i\right)\left[\psi^1(\alpha+\sum_j x_j)-\psi^1(\alpha+\beta+n)\right]+\left(n-\sum_i x_i\right)\left[\psi^1(\beta+n-\sum_j x_j)-\psi^1(\alpha+\beta+n)\right] \end{align*}\]

ただし，\(\psi^1(z)\)は1階のポリガンマ関数で，\(\Gamma(z)\)がガンマ関数のとき， \[\psi^1(z)=\frac{d^2}{dz^2} \log\Gamma(z)\]

GFV<-function(data) {
  (sum(data))*(psigamma(a+sum(data),1)-psigamma(a+b+length(data),1))+
    (length(data)-sum(data))*(psigamma(b+length(data)-sum(data),1)-psigamma(a+b+length(data),1))
}
GFV(ex)/length(ex)
## [1] 0.03259217

「WAIC(Watanabe-Akaike Widely Applicable Information Criterion)」は，経験損失と汎関数分散を\(n\)で割った値の和として定義される．

\[\mathrm{WAIC}=T_n+\frac{V_n}{n}=-\frac{1}{n}\sum_i\left\{\log E_\theta[p(x_i|\theta)]+E_\theta[(\log p(x_i|\theta))^2]-E_\theta[\log p(x_i|\theta)]^2\right\}\]

WAIC<-function(data) {
  T_loss(data)+GFV(data)/length(data)
}
WAIC(ex)
## [1] 0.6898985

理論的に，次のことが成り立つ（渡辺 2012: 118）． \[E_{X^n}(\mathrm{WAIC})=E_{X^n}(G_n)+o(1/n)\] このことは，実現可能性，正則条件を満たすか満たさないかにかかわらず，より一般的に成立する．スゲェ！

StanでWAIC

確認のために，stanでも具体例モデルを推定し，looパッケージでWAICを求めよう．

data {
  int N;
  int X[N];
}

parameters {
  real<lower=0,upper=1> q;
}

model {
  for (n in 1:N) {
    X[n] ~ bernoulli(q);
  }
}

generated quantities {
  real log_lik[N];
  for (n in 1:N) {
    log_lik[n] = bernoulli_lpmf(X[n]|q);
  }
}

library(rstan)
library(loo)
fit <- sampling(model1, data=list(X=ex,N=length(ex)),
                refresh=0,show_messages=FALSE,seed=1234)
## 実際には保存しているRDSを読み込み
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Warning: package 'loo' was built under R version 3.5.1
## This is loo version 2.0.0.
## **NOTE: As of version 2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. Visit mc-stan.org/loo/news for details on other changes.
stan_trace(fit,"q")

stan_hist(fit,"q")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

waic_est<-waic(extract(fit)$log_lik)$estimates
waic_est
##              Estimate         SE
## elpd_waic -20.6852365 1.48354512
## p_waic      0.9673747 0.09666489
## waic       41.3704729 2.96709024

looにおけるp_waicが汎関数分散に対応する．

c(waic_est[2,1],GFV(ex))
## [1] 0.9673747 0.9777652

また，looにおけるWAICは\(2n\)をかけた値になっている(Vehtari and Gelman 2014)．具体例データにおけるstanのWAIC，定義によるWAIC，汎化損失を比較する．

c(waic_est[3,1]/(2*n),WAIC(ex),gen_loss(ex))
## [1] 0.6895079 0.6898985 0.6232513

最後に，WAICが平均的に汎化損失のよい推定量になっているかを確認しよう．

## mean stan WAIC
stan_waic<-function(data) {
  fit <- sampling(tempmodel, data=list(X=data,N=length(data)),
              seed=1234,refresh=0,show_messages=FALSE)
  waic_est<-waic(extract(fit)$log_lik)$estimates
  waic_est[3,1]/(2*length(data))
}
nrep<-10000
stan_waic_dist<-numeric(nrep)
for (i in 1:nrep) {
  stan_waic_dist[i]<-stan_waic(rbinom(n,1,q))
} ##実際には保存しているRDSを読み込み
## Mean WAIC
WAIC_dist<-replicate(10000,WAIC(rbinom(n,1,q)))
## Mean general loss
gen_loss_dist<-replicate(10000,gen_loss(rbinom(n,1,q)))
## compare
boxplot(stan_waic_dist,WAIC_dist,gen_loss_dist,
        names=c("stan WAIC/2n","WAIC","G_n"),
        col=c("#FF00007F","#0000FF7F","#00FF007F"))

c(mean(stan_waic_dist),mean(WAIC_dist),mean(gen_loss_dist))
## [1] 0.6264795 0.6261056 0.6262251

一致かくにん！　よかった．

(3) BIC, WBIC

周辺尤度

事後分布の分母部分は，パラメータとデータの同時確率\(p(X^n,\theta)=p(X^n|\theta)\varphi(\theta)\)を\(\theta\)について周辺化したものであって，事前分布と確率モデルにおける\(X^n\)の確率分布である．これを「周辺尤度(marginal likelihood)」という．

\[p(X^n)=\int p(X^n|\theta)\varphi(\theta)d\theta\]

ベルヌイ確率モデルの具体例においては，以下のようになる． \[p(x^n)=\frac{B(\alpha +\sum_i x_i,\beta+n-\sum_i x_i)}{B(\alpha,\beta)}\]

marginal_likelihood<-function(data){
  beta(a+sum(data),b+length(data)-sum(data))/beta(a,b)
}
marginal_likelihood(ex)
## [1] 5.905118e-10

自由エネルギー

周辺尤度の対数の符号を逆転させたものを（逆温度1の）「自由エネルギー(free energy)」といって\(F_n\)で表す．

\[F_n=-\log p(X^n)=-\log\int p(X^n|\theta)\varphi(\theta)d\theta\]

自由エネルギー\(F_n\)の\(X^n\)についての平均は，真の分布のエントロピーの\(n\)倍と真の分布と周辺尤度の（\(X^n\)についての）カルバック=ライブラ(KL)距離に分解できる（渡辺 2012: 8-9）． \[E_{X^n}(F_n)=-n\int \log q(x) q(x)dx+\int \log \frac{q(x^n)}{p(x^n)} q(x^n)dx^n\] 真の分布のエントロピーは変わらないとすると，真の分布と周辺尤度の平均的な遠さを表すKL距離が小さくなるほど，自由エネルギーの平均は小さくなる．そこで，自由エネルギーが小さくなるほど周辺尤度がよく真の分布を近似していると言えるだろう．

f_energy<-function(data) {
  -log(marginal_likelihood(data))
}
f_energy(ex)
## [1] 21.25003

BIC

真の分布が確率モデルに対して正則であり，事後分布が正規分布に近似できる場合，積分のラプラス近似の手法を用いて，自由エネルギー（対数周辺尤度）を近似的に求められる．それが，「BIC(Bayesian Information Criterion)」であり，最尤推定から得られる最大対数尤度を用いた単純な式として定義される．

\[\mathrm{BIC}=-\sum_i \log p(x_i|\hat{\theta}_{ML})+\frac{d}{2}\log n\] 具体例におけるBICは，以下のようになる． \[\mathrm{BIC}=-\sum_i \left[x_i \log \bar{x} +(n-x_i)\log(1-\bar{x})\right]+\frac{\log n}{2}\]

BIC<-function(data) {
  -maximum_log_likelihood(data)+(1/2)*log(length(data))
}
BIC(ex)
## [1] 21.41533

WBIC

正則でないモデルにおいても，自由エネルギーの近似値を導出できるのが「WBIC(Watanabe Widely Applicable Bayesian Information Criterion)」である(Watanabe 2013，「広く使えるベイズ情報量規準(WBIC)」)．これは，逆温度が\(\beta=1/\log n\)のときの事後分布によって，対数損失 \[L_n(\theta)=-\frac{1}{n}\sum_i \log p(x_i|\theta)\] の\(n\)倍について期待値をとったものとして定義される． \[\mathrm{WBIC}=\frac{\int nL_n(\theta)\prod_i p(x_i|\theta)^\beta \varphi(\theta)d\theta}{\int \prod_i p(x_i|\theta)^\beta \varphi(\theta)d\theta}\] このとき自由エネルギー\(F_n\)と\(\mathrm{WBIC}\)は\(\log n\)のオーダーまで同じ漸近挙動をもつ．

具体例におけるWBICは，以下のようになる（Mathematicaさん）． \[\mathrm{WBIC}=n\psi\left(\alpha+\beta+\frac{n}{\log n}\right)-\left(n-\sum_i x_i\right)\psi\left(\beta+\frac{n-\sum_i x_i}{\log n}\right)-\left(\sum_i x_i\right)\psi\left(\alpha+\frac{\sum_i x_i}{\log n}\right)\] ただし，\(\psi(z)\)はディガンマ関数で，\(\Gamma(z)\)がガンマ関数のとき， \[\psi(z)=\frac{d}{dz} \log\Gamma(z)\]

WBIC<-function(data) {
  n<-length(data)
  x<-sum(data)
  n*digamma(a+b+n/log(n))-(n-x)*digamma(b+(n-x)/log(n))-x*digamma(a+x/log(n))
}
WBIC(ex)
## [1] 21.17431

最後に，BIC，WBIC，自由エネルギーの平均を比較する．

## Mean BIC
BIC_dist<-replicate(10000,BIC(rbinom(n,1,q)))
## Mean WBIC
WBIC_dist<-replicate(10000,WBIC(rbinom(n,1,q)))
## Mean free energy
free_energy_dist<-replicate(10000,f_energy(rbinom(n,1,q)))
## compare
boxplot(BIC_dist,WBIC_dist,free_energy_dist,
        names=c("BIC","WBIC","F_n"),
        col=c("#FF00007F","#0000FF7F","#00FF007F"))

c(mean(BIC_dist),mean(WBIC_dist),mean(free_energy_dist))
## [1] 19.53936 19.33736 19.43243

やったね！

文献