Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる（第3章）

ふと、Kéry & Schaubの“Bayesian Population Analysis using WinBUGS”のBUGSコードをStanに移植してみようとおもったので やってみました。まずは3章から。

なお、この本のBUGSコードとデータは、同書のサポートページComplete code and data files of the bookにて公開されています。

第3章は“Introduction to the Generalized Linear Model: The Simplest Model for Count Data”ということで、カウントデータのGLMです。まずは、3.3.1節の、仮想データをつかったポアソンGLMから。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## Specify model in Stan
GLM_Poisson_code <- "
data {
  int<lower=0> n;
  int<lower=0> C[n];
  real year[n];
}

parameters {
  real alpha;
  real beta1;
  real beta2;
  real beta3;
}

transformed parameters {
  vector[n] log_lambda;

  for (i in 1:n)
    log_lambda[i] <- alpha + beta1 * year[i] +
                     beta2 * year[i]^2 + beta3 * year[i]^3;
}

model {
  // Priors
  alpha ~ uniform(-20, 20);
  beta1 ~ uniform(-10, 10);
  beta2 ~ uniform(-10, 10);
  beta3 ~ uniform(-10, 10);

  // Likelihood
  C ~ poisson_log(log_lambda);
}

generated quantities {
  vector[n] lambda;

  lambda <- exp(log_lambda);
}
"

## Generate simulated data
set.seed(123)
data <- data.fn()

## Bundle data
stan.data <- list(C = data$C, n = length(data$C), year = data$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")

## 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)

data.fn()は、シミュレーションで仮想データを生成する関数で、同書のRコードそのままです。

このデータは、WinBUGSでは即エラーとなるのですが（year[i]^3あたりの数値が大きくなりすぎる）、Stanでも収束がよくありません。yearを標準化したデータをあたえたのが以下になります。

## 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)

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

こちらはすばやく収束しますので、Stanでも変数のスケーリングは大事ということがわかりました。

つづいて、3.5.1節の二項GLMです。

## Specify model in Stan
GLM_Binomial_code <- "
data {
  int<lower=0> nyears;
  int<lower=0> C[nyears];
  int<lower=0> N[nyears];
  real year[nyears];
}

parameters {
  real alpha;
  real beta1;
  real beta2;
}

transformed parameters {
  real logit_p[nyears];

  for (i in 1:nyears)
    // link function and linear predictor
    logit_p[i] <- alpha + beta1 * year[i] +
                  beta2 * pow(year[i],2);
}

model {
  // Priors
  alpha ~ normal(0, 100);
  beta1 ~ normal(0, 100);
  beta2 ~ normal(0, 100);

  // Likelihood
  for (i in 1:nyears){
    // 1. Distribution for random part
    C[i] ~ binomial_logit(N[i], logit_p[i]);
   }
}

generated quantities {
  real<lower=0,upper=1> p[nyears];

  for (i in 1:nyears)
    p[i] <- inv_logit(logit_p[i]);
}
"

## Generate simulated data
set.seed(123)
data <- data.fn(nyears = 40, alpha = 1, beta1 = -0.03, beta2 = -0.9)

## Bundle data
stan_data <- list(C = data$C, N = data$N,
                  nyears = length(data$C), year = data$YR)

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

## Parameters monitored
params <- c("alpha", "beta1", "beta2", "p")

## MCMC settings
ni <- 2500
nt <- 2
nb <- 500
nc <- 3

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

こちらも問題なく収束します。

このあたりはまだ離散値パラメーターもないので、すんなりと移植できました。

タグ：STAn