階層モデルと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
			

				

					

				

					
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
				
						
## 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]);
  }
}
			

				

					

				

					
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
				
						
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);
  }
}
			

				

					

				

					
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
				
						
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);
  }
}