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

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

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

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

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

１．モデル評価には２種類ある。AICとBICは見てるところが違うよ。

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

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

の３点です。

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

２年ぐらい前に、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);
}

data{

int N;

int Y;

real alpha;

real beta;

}

parameters{

real 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);
  }
}

data{

int N;

int T;

int Y[N];

}

parameters{

real 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]);
  }
}

data{

int N;

int T;

int Y[N];

}

parameters{

vector[N] theta;

real mu;

real 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);
  }
}

data{

int N;

int T;

int Y[N];

}

parameters{

real mu;

real 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)

100

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)