Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる（第7章その2）：Taglibro de H：SSブログ

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

Bayesian Population ..｜Bayesian Population .. ブログトップ

Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる（第7章その2）　[統計] [編集]

第7章その1のつづきです。

7.6. Models with time and group effects / 7.6.1. Fixed group and time effects

cjs_add_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
  int<lower=1> g;
  int<lower=1,upper=g> group[nind];
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];
  real beta1;

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
  beta1 <- 0;      // Corner constraint
}

parameters {
  real<lower=0,upper=1> mean_phi;
  real<lower=0,upper=1> mean_p;
  real gamma[n_occasions - 1];
  real<lower=0,upper=1> p_g[g];
  real beta2;
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;
  vector[g] beta;

  beta[1] <- beta1;
  beta[2] <- beta2;

  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- inv_logit(beta[group[i]] + gamma[t]);
      p[i, t] <- p_g[group[i]];
    }
  }

  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  mean_phi ~ uniform(0, 1);    // Prior for mean survival
  mean_p ~ uniform(0, 1);      // Prior for mean recapture
  beta2 ~ normal(0, 10)T[-10,10]; // Prior for difference in male
                                  // and female survival
  p_g ~ uniform(0, 1);         // Priors for group-spec. recapture
  gamma ~ normal(0, 10);       // Priors for time effects

  // Likelihood
  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}

generated quantities {
  real<lower=0,upper=1> phi_g1[n_occasions - 1];
  real<lower=0,upper=1> phi_g2[n_occasions - 1];

  for (t in 1:(n_occasions - 1)) {
    phi_g1[t] <- inv_logit(gamma[t]); // Back-transformed survival of males
    phi_g2[t] <- inv_logit(gamma[t] + beta[2]); // Back-transformed survival of females
  }
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2],
                  g = length(unique(group)),
                  group = group)

## Initial values
inits <- list(list(gamma = rnorm(n.occasions - 1),
                   beta = c(0, rnorm(1)),
                   p_g = runif(length(unique(group)), 0, 1)),
              list(gamma = rnorm(n.occasions - 1),
                   beta = c(0, rnorm(1)),
                   p_g = runif(length(unique(group)), 0, 1)),
              list(gamma = rnorm(n.occasions - 1),
                   beta = c(0, rnorm(1)),
                   p_g = runif(length(unique(group)), 0, 1)),
              list(gamma = rnorm(n.occasions - 1),
                   beta = c(0, rnorm(1)),
                   p_g = runif(length(unique(group)), 0, 1)))

## Parameters monitored
params <- c("phi_g1", "phi_g2", "p_g", "beta")

## MCMC settings
ni <- 3000
nt <- 2
nb <- 1000
nc <- 4

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

7.6.2. Fixed group and random time effects

群間誤差を多変量正規分布でモデル化しています。

cjs_temp_corr_code <- "
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }

  int last_capture(int[] y_i) {
    for (k_rev in 0:(size(y_i) - 1)) {
      int k;
      k <- size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }

  matrix prob_uncaptured(int nind, int n_occasions,
                         matrix p, matrix phi) {
    matrix[nind, n_occasions] chi;

    for (i in 1:nind) {
      chi[i, n_occasions] <- 1.0;
      for (t in 1:(n_occasions - 1)) {
        int t_curr;
        int t_next;

        t_curr <- n_occasions - t;
        t_next <- t_curr + 1;
        chi[i, t_curr] <- (1 - phi[i, t_curr]) +
                          phi[i, t_curr] *
                          (1 - p[i, t_next - 1]) *
                          chi[i, t_next];
      }
    }
    return chi;
  }
}

data {
  int<lower=0> nind;
  int<lower=2> n_occasions;
  int<lower=0,upper=1> y[nind, n_occasions];
  int<lower=1> g;
  int<lower=1,upper=g> group[nind];
  real df;
  matrix[g,g] R;
}

transformed data {
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];

  for (i in 1:nind)
    first[i] <- first_capture(y[i]);
  for (i in 1:nind)
    last[i] <- last_capture(y[i]);
}

parameters {
  vector<lower=0,upper=1>[g] mean_phi;
  real<lower=0,upper=1> p_g[g];
  matrix[n_occasions-1, g] eta_phi;
  cov_matrix[g] Sigma;
}

transformed parameters {
  matrix<lower=0,upper=1>[nind, n_occasions - 1] phi;
  matrix<lower=0,upper=1>[nind, n_occasions - 1] p;
  matrix<lower=0,upper=1>[nind, n_occasions] chi;
  vector[g] mu_phi;

  for (u in 1:g)
    mu_phi[u] <- logit(mean_phi[u]);

  // Constraints
  for (i in 1:nind) {
    for (t in 1:(first[i] - 1)) {
      phi[i, t] <- 0;
      p[i, t] <- 0;
    }
    for (t in first[i]:(n_occasions - 1)) {
      phi[i, t] <- inv_logit(eta_phi[t, group[i]]);
      p[i, t] <- p_g[group[i]];
    }
  }
  chi <- prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  // for survival parameters
  Sigma ~ inv_wishart(df, R);

  for (t in 1:(n_occasions - 1))
    increment_log_prob(multi_normal_log(row(eta_phi, t), mu_phi, Sigma));
  mean_phi ~ uniform(0, 1);   // Priors on mean group-spec. survival

  // for recapture parameters
  p_g ~ uniform(0, 1);        // Priors for group-spec. recapture

  // Likelihood
  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        1 ~ bernoulli(phi[i, t - 1]);
        y[i, t] ~ bernoulli(p[i, t - 1]);
      }
      1 ~ bernoulli(chi[i, last[i]]);
    }
  }
}
"

## Bundle data
stan_data <- list(y = CH,
                  nind = dim(CH)[1],
                  n_occasions = dim(CH)[2],
                  g = length(unique(group)),
                  group = group,
                  R = matrix(c(5, 0, 0, 1), ncol = 2),
                  df = 3)

## Initial values
inits <- function(){list(p_g = runif(length(unique(group)), 0, 1),
                         Omega = matrix(c(1, 0, 0, 1), ncol = 2))}  

## Parameters monitored
params <- c("eta_phi", "p_g", "Sigma", "mean_phi")

# MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 3

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

タグ：STAn