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)) }
するとこのようなものが出来る。
空間相関を考慮しないバージョンはこんな感じ。
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]); }
参考
stanコードの記述にあたっては下記サイトが参考になったというかほとんどそのまんまです。