Stan: Proper CAR model (再修正版)

さらに前回の修正です。 γの範囲がただしく指定されていなかったのを修正しました（文献はしっかりよむように > 自分）。

C, M, gamma_max, gamma_minをRで計算して、データとしてRStanにわたすようにしました。

##
## Proper CAR model
##
load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/Y.RData"))
N_site <- length(Y)

## Adjacent Matrix
W <- matrix(0, nrow = N_site, ncol = N_site)
for (i in 1:N_site) {
  for (j in 1:N_site) {
    if (abs(i - j) == 1)
      W[i, j] <- 1
  }
}
rsum <- apply(W, 1, sum)
C <- W / rsum
M <- diag(1 / rsum)
M2 <- svd(M)$u %*% diag(sqrt(svd(M)$d)) %*% t(svd(M)$v)    # M^(1/2)
gamma_max <- 1 / max(eigen(solve(M2) %*% C %*% M2)$values)
gamma_min <- 1 / min(eigen(solve(M2) %*% C %*% M2)$values)

library(rstan)
library(parallel)

stan_code <- "
data {
  int<lower=0> N_site;
  int<lower=0> Y[N_site];
  matrix[N_site, N_site] C;
  matrix[N_site, N_site] M;
  real gamma_max;
  real gamma_min;
}
transformed data {
  vector[N_site] mu;
  matrix[N_site, N_site] I;

  for (i in 1:N_site) {
    mu[i] <- 0;
    for (j in 1:N_site) {
      I[i, j] <- if_else(i == j, 1, 0);
    }
  }
}
parameters {
  vector[N_site] r;
  real beta;
  real<lower=0> s2;
  real<lower=gamma_min,upper=gamma_max> gamma;
}
transformed parameters {
  cov_matrix[N_site] V;

  V <- s2 * inverse(I - gamma * C) * M;
}
model {
  Y ~ poisson(exp(beta + r));
  r ~ multi_normal(mu, V);
  beta ~ normal(0, 100);
  s2 ~ uniform(0, 100);
  gamma ~ uniform(gamma_min, gamma_max);
}
generated quantities {
  vector[N_site] Y_hat;

  Y_hat <- exp(beta + r);
}
"

seeds <- c(12, 123, 1234, 12345)
inits <- list()
inits[[1]] <- list(beta = 1.0, r = rep(0, N_site), s2 = 1.0, gamma = 0.1)
inits[[2]] <- list(beta = 1.5, r = rep(0, N_site), s2 = 0.5, gamma = 0.2)
inits[[3]] <- list(beta = 2.0, r = rep(0, N_site), s2 = 0.2, gamma = 0.3)
inits[[4]] <- list(beta = 2.5, r = rep(0, N_site), s2 = 0.1, gamma = 0.4)
model <- stan_model(model_code = stan_code)
fit.list <- mclapply(1:4,
                     function(i) {
                       sampling(model,
                                data = list(Y = Y, N_site = N_site,
                                            C = C, M = M,
                                            gamma_max = gamma_max,
                                            gamma_min = gamma_min),
                                pars = c("beta", "r", "s2", "gamma", "Y_hat"),
                                chains = 1,
                                iter = 6000, warmup = 1000, thin = 5,
                                seed = seeds[i], init = list(inits[[i]]))
                     }, mc.cores = getOption("mc.cores", 4L))
fit <- sflist2stanfit(fit.list)

print(fit)
traceplot(fit, pars = c("beta", "s2", "gamma"))

pars <- extract(fit)

## plot
library(ggplot2)
df <- data.frame(Site = 1:N_site,
                 Y = Y,
                 Y_mean = apply(pars2$Y_hat, 2, mean),
                 Y_90 = apply(pars2$Y_hat, 2, quantile, probs = 0.9),
                 Y_10 = apply(pars2$Y_hat, 2, quantile, probs = 0.1))
ggplot(df) +
  geom_point(aes(x = Site, y = Y)) +
  geom_line(aes(x = Site, y = Y_mean)) +
  geom_ribbon(aes(x = Site, ymin = Y_10, ymax = Y_90), alpha = 0.5)

タグ：STAn