読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

階層ベイズでもサンプルサイズを増やしたらベイズ信頼区間の幅は細くなってくれるのか

系列長 N の時系列データがあるとして, これに以下のようなモデルを当てはめます.

 \mu_i  \sim N(\mu_{i-1},\sigma^2_2)
 x_i  \sim N(\mu_i,\sigma^2_1)

階層事前分布として \sigma_1, \sigma_2 には平均 10 の指数分布, \mu_1 には幅の広い一様分布を仮定します.

いま \mu_i の95%ベイズ信頼区間を求めたいとします.

ベイズ信頼区間の幅はサンプルサイズ(系列長 N)を増やしたとき, せまくなってくれるでしょうか.

  • サンプルサイズが増えるのだから, せまくなる
  • パラメータ \mu_i はサンプルサイズと同じオーダーで増えるのだから, サンプルサイズを増やしてもせまくならない

どちらもありそうな気がします.

理論的に考えるのはむずかしそうなので, 数値的に実験してみます.

MCMC(Stan)を使ってパラメータの分布を推定し, 得られた MCMC サンプルの2.5%点と97.5%点を取ることで, ベイズ信頼区間を求めました. 下の図です.

f:id:abrahamcow:20161102023714p:plain

乱数で x_i を生成しています.

赤い線が最初の10個のサンプルのみを使った場合のベイズ信頼区間, 青い線が50個使ったとき, 灰色が100個すべてを使ったときです.

見ての通り, 信頼区間の幅はサンプルサイズが増えるごとにせまくなっていますね.

  • Q. 階層ベイズでもサンプルサイズを増やしたらベイズ信頼区間の幅は細くなるのか.
  • A. なる.

ついでに \sigma_1の事後分布は以下のように推定されました.

f:id:abrahamcow:20161102024351p:plain

先ほどと同様, 赤い線が最初の10個, 青い線が50個を使ったとき, 黒が100個すべてを使ったときです.

サンプルサイズが増えるごとに形が急峻になっていきます.

\sigma_2 の事後分布は以下のようになりました.

f:id:abrahamcow:20161102024817p:plain

シミュレーションに使用した R のコードは以下です.

library(rstan)
library(cowplot)
set.seed(1234)
x <-cumsum(rnorm(100))
plot(x,type="l")
rstan_options(auto_write = TRUE)
locallevelmodel <-stan_model("locallevelmodel.stan")

dat10 <- list(x=x[1:10],N=10)
fit10 <-sampling(locallevelmodel,dat10,chain=1,iter=10000,control=list(adapt_delta=0.95))
ex10 <-extract(fit10)
lower10=apply(ex10$mu, 2,quantile,0.025)
upper10=apply(ex10$mu, 2,quantile,0.975)

dat50 <- list(x=x[1:50],N=50)
fit50 <-sampling(locallevelmodel,dat50,chain=1,iter=10000,control=list(adapt_delta=0.9))
ex50 <-extract(fit50)
lower50=apply(ex50$mu, 2,quantile,0.025)
upper50=apply(ex50$mu, 2,quantile,0.975)

dat100 <- list(x=x,N=100)
fit100 <-sampling(locallevelmodel,dat100,chain=1,iter=10000)
ex100 <-extract(fit100)
lower100=apply(ex100$mu, 2,quantile,0.025)
upper100=apply(ex100$mu, 2,quantile,0.975)

ggplot()+
  geom_line(aes(x=1:100,y=x))+
  geom_ribbon(aes(x=1:100,ymin=lower100,ymax=upper100),fill=NA,colour="grey60")+
  geom_ribbon(aes(x=1:50,ymin=lower50,ymax=upper50),fill=NA,colour="blue")+
  geom_ribbon(aes(x=1:10,ymin=lower10,ymax=upper10),fill=NA,colour="red")+
  xlab("")

ggplot()+
  geom_density(aes(x=ex100$sigma1))+
  geom_density(aes(x=ex50$sigma1),colour="blue")+
  geom_density(aes(x=ex10$sigma1),colour="red")

ggplot()+
  geom_density(aes(x=ex100$sigma2))+
  geom_density(aes(x=ex50$sigma2),colour="blue")+
  geom_density(aes(x=ex10$sigma2),colour="red")+
  xlab("")

推定に用いた Stan のコードは以下です.

data{
  int N;
  real x[N];
}
parameters{
  real mu[N];
  real<lower=0> sigma1;
  real<lower=0> sigma2;
}
model{
  mu[2:N] ~ normal(mu[1:(N-1)],sigma2);
  x ~ normal(mu,sigma1);
  sigma1 ~ exponential(0.1);
  sigma2 ~ exponential(0.1);
}