階層ベイズと自由エネルギー

先日、広島ベイズ塾の春合宿がありました。

そこで発表した、「階層ベイズと自由エネルギー」の資料をアップしました。

自由エネルギーというのは、負の対数周辺尤度のことで、ベイズファクターの計算で使う周辺尤度を対数とって‐1をかけたものです。

スライド内容を要約すると、

1.モデル評価には2種類ある。AICとBICは見てるところが違うよ。

2.階層ベイズではWAICはどういう予測をするかで値が変わるが、自由エネルギーは変わらないよ。

3.心理学では自由エネルギーのほうが知りたい値かもしれないね。

の3点です。

 

 

2年ぐらい前に、HijiyamaRで階層ベイズとWAICについて発表したものがありますが、それの続きになります。

ただ、内容的には松浦健太郎さんのブログ記事、「階層ベイズモデルとWAIC」のほうが断然わかりやすいので、こちらを先に見てもらったほうがいいかもしれません。一応、上のスライドだけみてもわかるようにはなってると思います。

想定読者はこれから統計モデリングを使っていこうとしている心理学者ですが、心理学以外の人にもわかる内容ですのでよかったら見てみてください。

あと、数学的なところ、自信がないので間違えていたらご指摘いただけると助かります。

スライドは図を張り付けてしまってるので、コードのコピペができません。なので、コードだけを下に貼っています。こちらもあわせてどうぞ。

model0.stan

data{
  int N;
  int Y;
  real alpha;
  real beta;
}

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

model{
  target += binomial_lpmf(Y|N,theta);
  target += beta_lpdf(theta|alpha,beta);
}

 

model1.stan

data{
  int N;
  int T;
  int Y[N];
}

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

model{
  target += binomial_lpmf(Y|T,theta);
}

generated quantities{
  vector[N] log_lik;
  for(n in 1:N){
    log_lik[n] = binomial_lpmf(Y[n]|T,theta);
  }
}

 

model2.stan

data{
  int N;
  int T;
  int Y[N];
}

parameters{
  vector<lower=0,upper=1>[N] theta;
  real<lower=0,upper=1> mu;
  real<lower=0> sigma;
}

model{
  target += binomial_lpmf(Y|T,theta);
  target += beta_lpdf(theta|mu/sigma, (1-mu)/sigma);
  target += student_t_lpdf(sigma|4,0,1);
}

generated quantities{
  vector[N] log_lik;
  for(n in 1:N){
    log_lik[n] = binomial_lpmf(Y[n]|T,theta[n]);
  }
}

 

model3.stan

data{
  int N;
  int T;
  int Y[N];
}

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

model{
  target += beta_binomial_lpmf(Y|T,mu/sigma, (1-mu)/sigma);
  target += student_t_lpdf(sigma|4,0,2.5);
}

generated quantities{
  vector[N] log_lik;
  for(n in 1:N){
    log_lik[n] = beta_binomial_lpmf(Y[n]|T,mu/sigma, (1-mu)/sigma);
  }
}

Rコード

library(rstan)
library(bridgesampling)
library(loo)
library(gamlss)

##function
waic2 <- function(log_lik){
  lppd <- sum(log(colMeans(exp(log_lik))))
  p_waic <- sum(apply(log_lik,2,var))
  waic <- -2*lppd+2*p_waic
  return(list(waic=waic,lppd=lppd,p_waic=p_waic))
}

binomial <- function(theta){
  temp <- function(theta) dbinom(Y,N,theta)*dbeta(theta,alpha,beta)
  return(sapply(theta,temp))
}

##bridge-sampling test

#data
N <- 20
Y <- 12
alpha <- 1
beta <- 1

#integrate
log(integrate(binomial,0,1)$value)

#stan and bridge-sampling
model0 <- stan_model("model0.stan")
data <- list(N=N,Y=Y,alpha=alpha,beta=beta)
fit0 <- sampling(model0,data=data,iter=20000,chains=2,cores=2)
bs0 <- bridge_sampler(fit0)
logml(bs0)


##hierarchical model

#model
model1 <- stan_model("model1.stan")
model2 <- stan_model("model2.stan")
model3 <- stan_model("model3.stan")

#data
rstan::expose_stan_functions("stan_func.stan")
Y <- betabinom_rng(100,20,0.5,0.1)

rbetabinom <- function(N,Tr,mu,sigma){
  theta <- rbeta(N,mu/sigma,(1-mu)/sigma)
  return(rbinom(N,Tr,theta))
}
set.seed(123)
N <- 100
Tr <- 20
mu <- 0.5
sigma <- 0.1
Y <- rbetabinom(N,Tr,mu,sigma)
hist(Y)


data <- list(N=N,T=20,Y=Y)

#sampling
fit1 <- sampling(model1,data=data,iter=21000,warmup=1000,chains=2,cores=2)
fit2 <- sampling(model2,data=data,iter=21000,warmup=1000,chains=2,cores=2)
fit3 <- sampling(model3,data=data,iter=21000,warmup=1000,chains=2,cores=2)

print(fit1,pars="theta")
print(fit2,pars=c("mu","sigma"))
print(fit3,pars=c("mu","sigma"))

#WAIC
waic2(rstan::extract(fit1)$log_lik)
waic2(rstan::extract(fit2)$log_lik)
waic2(rstan::extract(fit3)$log_lik)

#marginal likelihood
bs1 <- bridge_sampler(fit1)
-2*logml(bs1)
bs2 <- bridge_sampler(fit2)
-2*logml(bs2)
bs3 <- bridge_sampler(fit3)
-2*logml(bs3)

exp(logml(bs2)-logml(bs1))

#binomial regression(ML)
result <- glm(cbind(Y,rep(Tr,N)-Y)~1,
              family=binomial(link="identity"))
AIC(result)

#beta-binomial regression(ML)
result <- gamlss(cbind(Y,rep(Tr,N)-Y)~1,
                 family=BB(mu.link = "identity"))
AIC(result)

BIC(result)

summary(result)

 

 

This entry was posted in 心理統計学, Stan. Bookmark the permalink.