「データ解析のための統計モデリング入門」第11章メモ

p.254の図11.7をstanを使って再現したかった。

stanコードはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
}

model {
   for (j in 2:N)
      r[j] ~ normal(r[j-1], s_r); // 1次のCAR model
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(r[Idx_obs[j]]);
}

generated quantities { 
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(r[j]);
}

R側のコードはこんな感じ。

# データを欠測させる
no <- c(6, 9, 12, 13, 26:30)
Idx_obs <- (1:50)[-no]
Y_obs <- Y[Idx_obs]
N <- length(Y)
N_obs <- length(Y_obs)

# サンプリング
dlist <- list(N = N, N_obs = N_obs, Idx_obs = Idx_obs, Y_obs = Y_obs)
fit2 <- stan("sec11_2.stan", data = dlist)

# 結果から中央値と信頼区間を算出
pred <- extract(fit2)
ci <- apply(pred$Y_mean, 2, quantile, probs = c(0.5, 0.025, 0.975))

# プロット
plot(Y, type = "n", ylim = c(0, 25))
polygon(c(1:N, N:1), 
        c(ci[2, ], rev(ci[3, ])),
        border = NA,
        col = "gray90")
points(Y)
points(no, Y[no], pch = 19)
lines(m, lty = 2)
lines(ci[1, ])
for(n in no){
  polygon(c(n-0.5, n+0.5, n+0.5, n-0.5),
          c(0, 0, 100, 100), 
          border = NA, 
          col = rgb(.3, .3, .3, .3))
}

するとこのようなものが出来る。
f:id:Rion778:20161008220158p:plain
空間相関を考慮しないバージョンはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
   real mu;
}

model {
   for (j in 1:N)
      r[j] ~ normal(0, s_r);
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(mu + r[Idx_obs[j]]);
}

generated quantities {
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(mu + r[j]);
}

f:id:Rion778:20161008220315p:plain

参考

stanコードの記述にあたっては下記サイトが参考になったというかほとんどそのまんまです。