階層モデルとWBIC

階層モデルだとWBICの推定値が上手くいかないようなのですが、
周辺化する以外にいい方法を知っている人がいたら教えてもらえると嬉しいです。

またシミュレーションを追加して資料を修正するかもしれません。

コード

上のスライドで使ったRコードは以下です。

## library ---------------------
library(rstan)
library(magrittr)
library(bridgesampling)

## functions -------------------

wbic <- function(log_lik){
  -mean(rowSums(log_lik))
}

wbic_v_t <- function(log_lik){
  gf_variance <- function(log_lik){
    sum(colMeans(log_lik^2)-colMeans(log_lik)^2)
  }
  N <- ncol(log_lik)
  return(wbic(log_lik)+(1/(2*log(N)))*gf_variance(log_lik))
}

Fn.nb <- function(x,a,b){
  temp <- function(theta){
    lik <- 1
    for(i in 1:length(x)){
      lik <- lik * dnbinom(x[i],size=a,prob=(1/(1+theta)))
    }
    lik*dexp(theta,b)
  }
  integrate(temp,lower=0,upper=Inf)$value %>% log %>% multiply_by(-1)
}


## dataset ---------------------

N <- 150
set.seed(57)
a <- 3
b <- 1
theta_prime <- 1
Y <- rnbinom(N,size = a, prob = 1/(1+theta_prime))
Y %>% hist

standata <- list(N=N,Y=Y,a=a, b=b)



## stan setting --------------------------------

model.p <- stan_model("poisson_h.stan")
model.p.WBIC <- stan_model("poisson_h_WBIC.stan")
model.nb <- stan_model("negbinom.stan")
model.nb.WBIC <- stan_model("negbinom_WBIC.stan")
model.int.WBIC <- stan_model("poisson_h_int_WBIC.stan")


## Free energy --------------------------

Fn.nb(dat$Y,a,b)

## poisson hierarchical ------------------------

# bridesampling
fit.p <- sampling(model.p,
                data=standata,
                iter=11000,
                warmup=1000,
                chains=4,
                cores=4)

print(fit.p,pars=c("theta"))

bs_obj <- bridge_sampler(fit.p)
bs_obj %>% logml %>% multiply_by(-1)


#WBIC
fit.p.WBIC <- sampling(model.p.WBIC,
                     data=standata,
                     iter=3000,
                     warmup=1000,
                     chains=4,
                     cores=4)

fit.p.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic

fit.p.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic_v_t


## negative binomail (non hierarchical) --------------------

#bridesampling
fit.nb <- sampling(model.nb,
                   data=standata,
                   iter=11000,
                   warmup=1000,
                   chains=4,
                   cores=4,
                   show_messages = FALSE)


print(fit.nb,pars=c("theta"))

bs_obj.nb <- bridge_sampler(fit.nb)
bs_obj.nb %>% logml %>% multiply_by(-1)

#WBIC
fit.nb.WBIC <- sampling(model.nb.WBIC,
                        data=standata,
                        iter=3000,
                        warmup=1000,
                        chains=4,
                        cores=4,
                        verbose = TRUE)

fit.nb.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic

fit.nb.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic_v_t


## integraged model ------------------------------

# WBIC
fit.int.WBIC <- sampling(model.int.WBIC,
                         data=standata,
                         iter=3000,
                         warmup=1000,
                         chains=4,
                         cores=4,
                         verbose = TRUE)

fit.int.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic

fit.int.WBIC %>% 
  rstan::extract() %$% 
  log_lik %>% 
  wbic_v_t

次に、Stanコードです。

bridgesamplingのため、すべてtarget記法で書いています。なお以下のコードはWBIC用のコードなので、bridgesamplingをする場合はmodelの1.0/log(N)をかけている部分を省いてください。

最初に階層モデルです。

data{
  int N;
  int Y[N];
  real a;
  real b;
}

parameters{
  real<lower=0> lambda[N];
  real<lower=0> theta;
}

model{
  target += 1.0/log(N) * poisson_lpmf(Y | lambda);
  target += gamma_lpdf(lambda | a, theta);
  target += exponential_lpdf(theta | b);
}

generated quantities{
  real log_lik[N];
  for(n in 1:N){
    log_lik[n] = poisson_lpmf(Y[n] | lambda[n]);
  }
}

次に、周辺化モデルです。

data{
  int N;
  int Y[N];
  real a;
  real b;
}

parameters{
  real<lower=0> theta;
}

model{
  target += 1/log(N) * neg_binomial_lpmf(Y | a, theta);
  target += exponential_lpdf(theta | b);
}

generated quantities{
  real log_lik[N];
  for(n in 1:N){
    log_lik[n] = neg_binomial_lpmf(Y[n] | a, theta);
  }
}