階層モデルだとWBICの推定値が上手くいかないようなのですが、
周辺化する以外にいい方法を知っている人がいたら教えてもらえると嬉しいです。
またシミュレーションを追加して資料を修正するかもしれません。
コード
上のスライドで使ったRコードは以下です。
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150 ## 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 <- 1for(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 <- 150set.seed(57)a <- 3b <- 1theta_prime <- 1Y <- rnbinom(N,size = a, prob = 1/(1+theta_prime))Y %>% histstandata <- 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 ------------------------# bridesamplingfit.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)#WBICfit.p.WBIC <- sampling(model.p.WBIC,data=standata,iter=3000,warmup=1000,chains=4,cores=4)fit.p.WBIC %>%rstan::extract() %$%log_lik %>%wbicfit.p.WBIC %>%rstan::extract() %$%log_lik %>%wbic_v_t## negative binomail (non hierarchical) --------------------#bridesamplingfit.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)#WBICfit.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 %>%wbicfit.nb.WBIC %>%rstan::extract() %$%log_lik %>%wbic_v_t## integraged model ------------------------------# WBICfit.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 %>%wbicfit.int.WBIC %>%rstan::extract() %$%log_lik %>%wbic_v_t
次に、Stanコードです。
bridgesamplingのため、すべてtarget記法で書いています。なお以下のコードはWBIC用のコードなので、bridgesamplingをする場合はmodelの1.0/log(N)をかけている部分を省いてください。
最初に階層モデルです。
123456789101112131415161718192021222324 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]);}}
次に、周辺化モデルです。
12345678910111213141516171819202122 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);}}