階層ベイズモデルのWAICを実感する
俺様用自習ノート
石田淳（関西学院大学社会学部）
2019年1月31日
要rstan,loo
はじめに
本ノートは，ベルヌイ確率モデルと，ベルヌイ分布のパラメータを生成するベータ分布，ベータ分布のパラメータの事前分布を仮定した階層モデルをもとに，WAICの意味を実感することを目的とする．
参考にしたのは，
- 渡辺澄夫「階層ベイズ法とWAIC」(http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/h_b.html)
- 松浦健太郎「階層ベイズモデルとWAIC」(http://statmodeling.hatenablog.com/entry/waic-with-hierarchical-models)
- 清水裕士「階層ベイズのWAICについて」(https://www.slideshare.net/simizu706/waic)
a<-2
b<-2
m<-10
n<-30ベータ分布からベルヌイ分布のパラメータが\(m\)個生成され， \[q_j\sim Beta(2,2)\]
それぞれの\(q_j\)から\(n\)個ずつベルヌイ確率変数\((X_{1j},\cdots,X_{nj})\)がサンプルされるケースを考える．
\[X_{ij} \sim Bern(q_j)\]
# example data
# q<-rbeta(m,a,b)
q<-c(0.8372668,0.8814993,0.2483013,0.1362599,0.3794142,
0.6132820,0.9144260,0.7322638,0.1465232,0.7190158)
# ex<-as.vector(sapply(q,function(x) rbinom(n,1,x)))
# write.csv(ex,"WAIC_H_model.csv")
ex<-read.csv("WAIC_H_model.csv")[,2]
ex## [1] 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 0 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 0 0 1 0 0 0 0 0 0 0
## [71] 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## [106] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 1 0 1 1 0
## [141] 1 0 0 0 0 0 0 1 0 0 0 1 0 1 1 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 0
## [176] 1 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1
## [211] 0 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0
## [246] 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 1 0
## [281] 0 0 1 1 0 0 0 1 1 1 1 0 1 0 1 1 1 1 1 1group<-as.vector(sapply(1:m,function(x) rep(x,n)))
ex_grouped<-sapply(1:m,function(x) sum(ex[((x-1)*30+1):(x*30)]))このモデルは，\(X_j=\sum_{i=1}^n X_{ij}\)のとき， \[X_j \sim Binom(n,q_j),\ \ \ \ \ \ q_j\sim Beta(2,2)\] あるいは，ベータ二項分布を用いて \[X_j \sim BetaBinom(n,2,2)\] と記述するモデルと（定数項\(\binom n x\)を除いて）一致する．
(1) 非階層モデルの汎化損失
事前分布\(q\sim Beta(1,1)\)，確率モデル\(X\sim Bern(q)\)を想定した非階層モデルを分析する．事後分布は， \[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^{\sum_i x_i}(1-q)^{n-\sum_i x_i}}{B(1+\sum_i x_i,1+n-\sum_i x_i)} \] である．予測分布は， \[p^*(X)=\int_0^1 p(X|q)p(q|x^n)dq=\left(\frac{1+\sum_i x_i}{1+1+n}\right)^X\left(\frac{1+n-\sum_i x_i}{1+1+n}\right)^{1-X}\] また，\(X\)の真の確率分布は， \[ \begin{aligned} q(X)&=\int_0^1 p(X|q)p(q)dq \\ &=\int_0^1 p(X,q)dq \\ &=\int_0^1 p(X)p(q)dq\ \ \ \mathrm{(by\ independency)} \\ &=\int_0^1 q\frac{q(1-q)}{B(2,2)}dq =\frac{1}{2} \end{aligned} \]
である．よって汎化損失は \[ \begin{aligned} G_n&=-E_X[\log p^*(X)] \\ &=-\frac{1}{2}\log\left(\frac{1+\sum_i x_i}{1+1+n}\right)-\frac{1}{2}\log\left(\frac{1+n-\sum_i x_i}{1+1+n}\right) \end{aligned} \] のはず・・・
gen_loss<-function(data,q) {
-q*log((1+sum(data))/(1+1+length(data)))-
(1-q)*log((1+length(data)-sum(data))/(1+1+length(data)))
}
gen_loss(ex,1/2)## [1] 0.6934982(2) 非階層モデルのWAIC
事前分布\(q\sim Beta(1,1)\)，確率モデル\(X\sim Bern(q)\)を想定した非階層モデルをstanで分析する．
data {
int N;
int X[N];
}
parameters {
real<lower=0,upper=1> q;
}
model {
for (n in 1:N) {
X[n] ~ bernoulli(q);
}
q ~ beta(1,1);
}
generated quantities {
real log_lik[N];
for (n in 1:N) {
log_lik[n] = bernoulli_lpmf(X[n]|q);
}
}
fit1 <- sampling(tempmodel1, data=list(X=ex,N=n*m),
seed=1234)##
## SAMPLING FOR MODEL '35558cde132dd64523e10f7aa6e8939b' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.088 seconds (Warm-up)
## Chain 1: 0.081 seconds (Sampling)
## Chain 1: 0.169 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '35558cde132dd64523e10f7aa6e8939b' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.086 seconds (Warm-up)
## Chain 2: 0.081 seconds (Sampling)
## Chain 2: 0.167 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '35558cde132dd64523e10f7aa6e8939b' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.086 seconds (Warm-up)
## Chain 3: 0.092 seconds (Sampling)
## Chain 3: 0.178 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '35558cde132dd64523e10f7aa6e8939b' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.084 seconds (Warm-up)
## Chain 4: 0.088 seconds (Sampling)
## Chain 4: 0.172 seconds (Total)
## Chain 4:stan_trace(fit1,"q")
stan_hist(fit1,"q")
WAIC（関数は松浦さんのウェブページのもの）を計算． \[ \mathrm{WAIC}=-\frac{1}{n*m}\sum_i\left\{\log E_q[p(x_i|q)]+E_q[(\log p(x_i|q))^2]-E_q[\log p(x_i|q)]^2\right\} \]
# http://statmodeling.hatenablog.com/entry/waic-with-hierarchical-models
waic <- function(log_likelihood) {
training_error <- - mean(log(colMeans(exp(log_likelihood))))
functional_variance_div_N <- mean(colMeans(log_likelihood^2) - colMeans(log_likelihood)^2)
waic <- training_error + functional_variance_div_N
return(waic)
}
waic(extract(fit1)$log_lik)## [1] 0.6962547looパッケージでも計算する．
loo::waic(extract(fit1)$log_lik)$estimates## Estimate SE
## elpd_waic -208.876676 0.435418928
## p_waic 1.038752 0.002839586
## waic 417.753352 0.870837855looのwaicは\(2mn\)をかけたもの．
2*m*n*waic(extract(fit1)$log_lik)## [1] 417.7528(3) 非階層モデルのWAIC2
事前分布\(q\sim Beta(1,1)\)，\(X=\sum X_i\)として，確率モデル\(X\sim Binom(n,q)\)を想定した非階層モデルをstanで分析する．
data {
int N ;
int X ;
}
parameters {
real<lower=0,upper=1> q;
}
model {
X ~ binomial(N, q);
q ~ beta(1,1);
}
generated quantities {
real log_lik;
log_lik = binomial_lpmf(X|N,q);// - log(choose(N,X));
}
fit2 <- sampling(tempmodel2, data=list(N=n*m,X=sum(ex)),
seed=1234)##
## SAMPLING FOR MODEL '9d20e81b1b7f49825275cd3038f7a34f' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.01 seconds (Warm-up)
## Chain 1: 0.01 seconds (Sampling)
## Chain 1: 0.02 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '9d20e81b1b7f49825275cd3038f7a34f' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.01 seconds (Warm-up)
## Chain 2: 0.011 seconds (Sampling)
## Chain 2: 0.021 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '9d20e81b1b7f49825275cd3038f7a34f' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.01 seconds (Warm-up)
## Chain 3: 0.011 seconds (Sampling)
## Chain 3: 0.021 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '9d20e81b1b7f49825275cd3038f7a34f' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.01 seconds (Warm-up)
## Chain 4: 0.01 seconds (Sampling)
## Chain 4: 0.02 seconds (Total)
## Chain 4:ここで，算出した1つだけのlog_likは \[\mathrm{log\_lik}=\log \binom nx + x\log q+(n-x)\log(1-q)\] なので，\(\mathrm{log\_lik}-\log \binom nx\)で，ベルヌイ確率モデルの対数尤度の総和と等しくなるはずである（ただし，stan上のgenerated quantitiesでは定数項の引き算ができない）．そこで，この項を引いて\(m*n\)で割ったものでwaicを計算する．
waic(as.matrix(extract(fit2)$log_lik-log(choose(n*m,sum(ex)))))/(m*n)## [1] 0.6954716loo::waic(as.matrix(extract(fit2)$log_lik-log(choose(n*m,sum(ex)))))$estimate## Warning: 1 (100.0%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.## Estimate SE
## elpd_waic -208.6416082 NA
## p_waic 0.4583116 NA
## waic 417.2832164 NA2*waic(as.matrix(extract(fit2)$log_lik-log(choose(n*m,sum(ex)))))## [1] 417.283対数尤度の総和を取る分，汎関数分散は小さくなるので（？），looのp_waicはBernモデルよりも小さくなる．また，サンプル1あつかいなので，SEはNAになる．
(3) 階層モデルのWAIC: それぞれのグループに1つサンプルが追加される場合，Bernを仮定
\[ X_i\sim Bern(q_{j(i)}), \\ q_j\sim Beta(a,b), \\ a\sim U(0,\infty), b\sim U(0,\infty) \]
data {
int N;
int X[N];
int K ;
int<lower=1,upper=K> Z[N] ;
}
parameters {
real<lower=0,upper=1> q[K];
real<lower=0> a ;
real<lower=0> b ;
}
model {
for (n in 1:N) {
X[n] ~ bernoulli(q[Z[n]]);
}
for (k in 1:K) {
q[k] ~ beta(a,b);
}
}
generated quantities {
real log_lik[N];
for (n in 1:N) {
log_lik[n] = bernoulli_lpmf(X[n]|q[Z[n]]);
}
}
fit3 <- sampling(tempmodel3, data=list(X=ex,N=m*n,Z=group,K=m),
seed=1234)##
## SAMPLING FOR MODEL 'f86853d2bb36f3befbe72a018a5fb4e6' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.286 seconds (Warm-up)
## Chain 1: 0.258 seconds (Sampling)
## Chain 1: 0.544 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f86853d2bb36f3befbe72a018a5fb4e6' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.279 seconds (Warm-up)
## Chain 2: 0.35 seconds (Sampling)
## Chain 2: 0.629 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f86853d2bb36f3befbe72a018a5fb4e6' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.293 seconds (Warm-up)
## Chain 3: 0.226 seconds (Sampling)
## Chain 3: 0.519 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f86853d2bb36f3befbe72a018a5fb4e6' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.272 seconds (Warm-up)
## Chain 4: 0.248 seconds (Sampling)
## Chain 4: 0.52 seconds (Total)
## Chain 4:stan_trace(fit3,"q")
stan_hist(fit3,"q")
stan_hist(fit3,pars=c("a","b"))
グループごとにサンプルが1つ追加されたときのWAIC．
\[ \mathrm{WAIC}_j=-\frac{1}{n_j}\sum_i\left\{\log E_{q_j}[p(x_{ij}|q_j)]+E_{q_j}[(\log p(x_{ij}|q_j))^2]-E_{q_j}[\log p(x_{ij}|q_j)]^2\right\} \]
g_waic<-sapply(1:10,function(x) waic(extract(fit3)$log_lik[,((x-1)*30+1):(x*30)]))
g_waic## [1] 0.4868960 0.4281674 0.3587832 0.1846904 0.6442055 0.6143413 0.4850669
## [8] 0.6707935 0.4252320 0.6897544それぞれのグループに同時に1つずつサンプルが追加されたときの予測の誤差の総和．
sum(g_waic)## [1] 4.987931looパッケージを使う．
loo::waic(extract(fit3)$log_lik)$estimates## Warning: 1 (0.3%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.## Estimate SE
## elpd_waic -149.640237 9.2541838
## p_waic 9.261679 0.9961579
## waic 299.280474 18.5083677looパッケージのwaicは，\(\sum_j 2n_j\mathrm{WAIC}_j\)．
sum(2*n*g_waic)## [1] 299.2758(4) 階層モデルのWAIC: それぞれのグループに1つサンプルが追加される場合，Binomを仮定
\[ X_j \sim Binom(n,q_j), \\ q_j\sim Beta(a,b), \\ a\sim U(0,\infty), b\sim U(0,\infty) \]
data {
int N ;
int M ;
int X[M] ;
}
parameters {
real<lower=0,upper=1> q[M];
real<lower=0> a ;
real<lower=0> b ;
}
model {
for (m in 1:M) {
X[m] ~ binomial(N, q[m]);
}
for (m in 1:M) {
q[m] ~ beta(a,b);
}
}
generated quantities {
real log_lik[M];
for (m in 1:M) {
log_lik[m] = binomial_lpmf(X[m]|N,q[m]) - log(choose(N,X[m]));
}
}
generated quantities内で\(\log \binom{n} {x_j}\)を引いている．
fit4 <- sampling(tempmodel4, data=list(N=n,M=m,X=ex_grouped),
seed=1234)##
## SAMPLING FOR MODEL 'f7fac907934be086bca0e499029cb846' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.109 seconds (Warm-up)
## Chain 1: 0.089 seconds (Sampling)
## Chain 1: 0.198 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'f7fac907934be086bca0e499029cb846' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.103 seconds (Warm-up)
## Chain 2: 0.109 seconds (Sampling)
## Chain 2: 0.212 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'f7fac907934be086bca0e499029cb846' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.11 seconds (Warm-up)
## Chain 3: 0.081 seconds (Sampling)
## Chain 3: 0.191 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'f7fac907934be086bca0e499029cb846' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.109 seconds (Warm-up)
## Chain 4: 0.129 seconds (Sampling)
## Chain 4: 0.238 seconds (Total)
## Chain 4:g_waic2<-sapply(1:m,function(x) waic(extract(fit4)$log_lik[,x,drop=F])/n)
g_waic2## [1] 0.4786904 0.4230205 0.3576310 0.1897819 0.6386872 0.6073514 0.4810702
## [8] 0.6626476 0.4263874 0.6837121sum(g_waic2)## [1] 4.94898loo::waic(extract(fit4)$log_lik)$estimates## Warning: 10 (100.0%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.## Estimate SE
## elpd_waic -148.470772 14.8108530
## p_waic 5.520904 0.3675923
## waic 296.941544 29.6217060sum(2*n*g_waic2)## [1] 296.9388(5) 階層モデルのWAIC: 新たなグループが追加される場合，BetaBinomを仮定
\(q_j\)を周辺化したベータ二項分布を仮定．
\[ X_j \sim BetaBinom(n,a,b), \\ a\sim U(0,\infty), b\sim U(0,\infty) \]
data {
int N ;
int M ;
int X[M] ;
}
parameters {
real<lower=0> a ;
real<lower=0> b ;
}
model {
for (m in 1:M) {
X[m] ~ beta_binomial(N, a, b);
}
}
generated quantities {
real log_lik[M];
for (m in 1:M) {
log_lik[m] = beta_binomial_lpmf(X[m]|N,a,b) - log(choose(N,X[m]));
}
}
generated quantities内で\(\log \binom{n} {x_j}\)を引いている．
fit5 <- sampling(tempmodel5, data=list(N=n,M=m,X=ex_grouped),
seed=1234)##
## SAMPLING FOR MODEL 'd30009487157dca0bdb8aa53b4650d86' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.118 seconds (Warm-up)
## Chain 1: 0.131 seconds (Sampling)
## Chain 1: 0.249 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'd30009487157dca0bdb8aa53b4650d86' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.115 seconds (Warm-up)
## Chain 2: 0.121 seconds (Sampling)
## Chain 2: 0.236 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'd30009487157dca0bdb8aa53b4650d86' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.122 seconds (Warm-up)
## Chain 3: 0.123 seconds (Sampling)
## Chain 3: 0.245 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'd30009487157dca0bdb8aa53b4650d86' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.116 seconds (Warm-up)
## Chain 4: 0.116 seconds (Sampling)
## Chain 4: 0.232 seconds (Total)
## Chain 4:stan_hist(fit3,pars=c("a","b"))
新たなグループが追加されたときのWAIC．
\[ \mathrm{WAIC}=-\frac{1}{m}\sum_j\left\{\log E_{(a,b)}[p(x_j|a,b)]+E_{(a,b)}[(\log p(x_j|a,b))^2]-E_{(a,b)}[\log p(x_j|a,b)]^2\right\} \]
waic(extract(fit5)$log_lik)## [1] 15.95994looパッケージを使う．
loo::waic(extract(fit5)$log_lik)$estimates## Warning: 1 (10.0%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.## Estimate SE
## elpd_waic -159.599844 13.5446276
## p_waic 1.906502 0.6514402
## waic 319.199688 27.0892552looパッケージのwaicは，\(2m\)倍したもの．
2*m*waic(extract(fit5)$log_lik)## [1] 319.1987まとめ
waic_tab<-cbind(loo::waic(extract(fit1)$log_lik)$estimates,
loo::waic(as.matrix(extract(fit2)$log_lik-log(choose(n*m,sum(ex)))))$estimate,
loo::waic(extract(fit3)$log_lik)$estimates,
loo::waic(extract(fit4)$log_lik)$estimates,
loo::waic(extract(fit5)$log_lik)$estimates)
colnames(waic_tab)<-c("non_H_bern","SE","non_H_binom","SE","H_grouped_bern","SE",
"H_grouped_binom","SE","H_new_group","SE")
waic_tab## non_H_bern SE non_H_binom SE H_grouped_bern
## elpd_waic -208.876676 0.435418928 -208.6416082 NA -149.640237
## p_waic 1.038752 0.002839586 0.4583116 NA 9.261679
## waic 417.753352 0.870837855 417.2832164 NA 299.280474
## SE H_grouped_binom SE H_new_group SE
## elpd_waic 9.2541838 -148.470772 14.8108530 -159.599844 13.5446276
## p_waic 0.9961579 5.520904 0.3675923 1.906502 0.6514402
## waic 18.5083677 296.941544 29.6217060 319.199688 27.0892552\(2\times サンプルサイズ\)でリスケールしたwaicであれば，階層モデルと非階層モデル，グループ内でのサンプル追加と，グループ追加のwaicを比較可能？ でもこれって結局何の指標？
属性でグループを構成する場合，各カテゴリがおおよそ可能な属性タイプを包含しているとき，新たな属性グループが出現することを想定するのは非現実的なので，各グループにそれぞれサンプルが追加される予測誤差を見ればよい？