R: Stanをためす (4) CARのような、でないような：Taglibro de H：SSブログ

	ブログをはじめるログイン

R: Stanをためす (3) 欠測値の..｜R: Stanをためす (5) 状態空間.. ブログトップ

R: Stanをためす (4) CARのような、でないような　[統計] [編集]

久保本11章のCARをやってみようとしたが、うまくいかなかった記録。

# KuboBook 11

library(coda)
library(parallel)

program <- "kubo11"

load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/Y.RData"))
N_site <- length(Y)
Adj <- c(sapply(2:(N_site - 1), function(a) c(a - 1, a + 1)))
Adj <- c(2, Adj, N_site - 1)
Weights <- rep(1, 2 * N_site - 2)
Num <- c(1, rep(2, N_site - 2), 1)
aindex1 <- rep(0, N_site)  # start
aindex2 <- rep(0, N_site)  # end
aindex1[1] <- 1
aindex2[1] <- Num[1]
for (i in 2:N_site) {
  aindex1[i] <- aindex2[i - 1] + 1
  aindex2[i] <- aindex2[i - 1] + Num[i]
}

# dump data
data.file <- paste(program, "_dump.R", sep = "")
dump(c("N_site", "Y", "Adj", "Weights", "aindex1", "aindex2"), data.file)

# initial values
init.file <- paste(program, "_init.R", sep = "")
beta <- 1
sigma <- 1
r <- rnorm(N_site, 0, 0.1)
dump(c("beta", "sigma", "r"), init.file)

# Stan
stan.home <- ("~/Documents/src/stan-src-1.0.2")
cur.dir <- getwd()
setwd(stan.home)
system(paste("make ", cur.dir, "/", program, sep = ""))
setwd(cur.dir)

warmup <- 100
iter <- 5000 + warmup
thin <- 5
seeds <- c(3141, 31415, 314159)
run <- function(chain) {
  cmd <- paste("./", program,
               " --data=", data.file,
               " --init=", init.file,
               " --seed=", seeds[chain],
               " --iter=", iter,
               " --warmup=", warmup, 
               " --thin=", thin,
               " --chain_id=", chain,
               " --samples=", program, ".chain", chain, ".csv", sep = "")
  system(cmd)
  s <- read.csv(paste(program, ".chain", chain, ".csv", sep = ""), comment.char = "#")
  mcmc(s, start = warmup + 1, thin = thin)
}

samples <- as.mcmc.list(mclapply(1:3, run, mc.cores = 4))

plot(samples[, c("beta", "sigma")])
summary(samples[, c("beta", "sigma", "r.1", "r.11")])

beta <- samples[, "beta"]
r <- lapply(1:N_site, function(i) samples[, paste("r", i, sep = ".")])

y <- sapply(1:N_site, function(i) unlist(beta) + unlist(r[i]))

## plot
plot(1:N_site, Y, type = "p", xlab = "Site", ylim = c(0, 25), las = 1)
lines(1:N_site, exp(apply(y, 2, mean)))
lines(1:N_site, exp(apply(y, 2, quantile, probs = 0.9)), lty = 2)
lines(1:N_site, exp(apply(y, 2, quantile, probs = 0.1)), lty = 2)

kubo11.stan

data {
  int<lower=0> N_site;
  int<lower=0> Y[N_site];
  int<lower=0> Adj[N_site * 2 - 2];
  real Weights[N_site * 2 - 2];
  int<lower=1> aindex1[N_site];
  int<lower=1> aindex2[N_site];
}
parameters {
  real r[N_site];
  real beta;
  real<lower=0> sigma;
}
transformed parameters {
  real s[N_site];
  real m[N_site];
  
  for (j in 1:N_site) {
    real c;
    real cr;

    c <- 0.0;
    cr <- 0.0;
    for (i in aindex1[j]:aindex2[j]) {
      c <- c + Weights[i];
      cr <- cr + Weights[i] * r[Adj[i]];
    }
    m[j] <- 0.8 * cr / c;
    s[j] <- sigma / sqrt(c);
  }
}
model {
  for (j in 1:N_site) {
    Y[j] ~ poisson(exp(beta + r[j]));
  }
  r ~ normal(m, s);
  beta ~ normal(0.0, 1.0e+2);
  sigma ~ uniform(0.0, 1.0e+4);
}

上の"0.8"が謎のパラメータ。しかしこれをいれないとrがどんどん発散してしまう。というあたりがうまくいかない。

グラフ