Posterior Predictive CheckをStanでやってみたメモ

Posterior Predictive Checkをやってみたメモです。Posterior Predictive Checkはモデルの適合度をしらべるもので、Bayesian Population Analysis using WinBUGS (BPA)の7.10節に説明があり、Bayesian Data Analysis in Ecology using Linear Models with R, BUGS, and Stanの10章、Bayesian Data Analysis, 3rd ed.の6.3節に解説があります。

Bayesian Population Analysis using WinBUGS: A hierarchical perspective

• 出版社/メーカー: Academic Press
• 発売日: 2011/10/11
• メディア: Kindle版

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

• 作者: Andrew Gelman
• 出版社/メーカー: Chapman and Hall/CRC
• 発売日: 2013/11/05
• メディア: ハードカバー

Bayesian Data Analysis in Ecology Using Linear Models with R, BUGS, and Stan

• 出版社/メーカー: Academic Press
• 発売日: 2015/04/04
• メディア: Kindle版

モデルとデータは、昨日のポアソンGLMのものをつかいました。Stanコードのうち、generated quantitiesブロックを以下のように修正します。

generated quantities {
  vector[n] lambda;
  int<lower=0> C_rep[n];
  real T_obs;
  real T_rep;

  lambda <- exp(log_lambda);
  T_obs <- 0;
  T_rep <- 0;
  for (i in 1:n) {
    C_rep[i] <- poisson_rng(lambda[i]);
    T_obs <- T_obs + (sqrt(C[i]) - sqrt(lambda[i]))^2;
    T_rep <- T_rep + (sqrt(C_rep[i]) - sqrt(lambda[i]))^2;
  }
}

MCMCの途中で計算される期待値を利用して、複製データセットを生成し、複製データセットと期待値との不一致測度を計算する。同様に、実際のデータと期待値との不一致測度も計算する。不一致測度には、モデルに応じてさまざまな統計量がつかわれる。ということで、このようなコードをかいたのですが、あまり自信はありません。不一致測度にはBPAの例題で使用されていたFreeman-Tukey統計量を使用しました。

以下のRコードでStanを実行しました。

## Bundle data
mean.year <- mean(data$year)             # Mean of year covariate
sd.year <- sd(data$year)                 # SD of year covariate

stan_data <- list(C = data$C, n = length(data$C),
             year = (data$year - mean.year) / sd.year)

## Initial values
inits <- function() list(alpha = runif(1, -2, 2),
                         beta1 = runif(1, -3, 3))

## Parameters monitored
params <- c("alpha", "beta1", "beta2", "beta3", "lambda",
            "T_obs", "T_rep")

## MCMC settings
ni <- 2000
nt <- 2
nb <- 1000
nc <- 3

## Stan
model <- stan_model(model_code = GLM_Poisson_code)
out <- sampling(model, data = stan_data, init = inits, pars = params,
                chains = nc, thin = nt, iter = ni, warmup = nb,
                seed = 1,
                open_progress = FALSE)

結果を図示します。

直線は1:1の線です。

> sum(samp$T_rep > samp$T_obs)/length(samp$T_obs)
[1] 0.7346667

理想的には、直線の上におちる点が50%になるということなのですが、どこかまちがっているかもしれません。