Bayesian Population Analysis using WinBUGSのコードをStanに移植してみる(第8章) [統計]
第7章がおわって、第8章 Estimation of Survival Using Mark-Recovery Data です。
8.2. The mark-recovery model as a state-space model
8.2.2. Analysis of a model with constant parameters
mr_ss_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; } real cell_prob(int n_occasions, row_vector s, row_vector r, int first, int last) { vector[n_occasions] pr; pr[first] <- (1 - s[first]) * r[first]; for (j in (first + 1):n_occasions - 1) pr[j] <- prod(segment(s, first, j - first)) * (1 - s[j]) * r[j]; for (j in 1:(first - 1)) pr[j] <- 0; for (t in 1:n_occasions - 1) pr[n_occasions] <- 1 - sum(head(pr, n_occasions - 1)); return pr[last]; } } data { int<lower=0> nind; int<lower=0> n_occasions; int<lower=0,upper=1> y[nind, n_occasions]; } 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 { real<lower=0,upper=1> mean_s; real<lower=0,upper=1> mean_r; } transformed parameters { matrix<lower=0,upper=1>[nind, n_occasions - 1] s; matrix<lower=0,upper=1>[nind, n_occasions - 1] r; // Constraints for (i in 1:nind) { for (t in 1:(first[i] - 1)) { s[i, t] <- 0; r[i, t] <- 0; } for (t in first[i]:n_occasions - 1) { s[i, t] <- mean_s; r[i, t] <- mean_r; } } } model { // Priors mean_s ~ uniform(0, 1); // Prior for mean survival mean_r ~ uniform(0, 1); // Prior for mean recapture // Likelihood for (i in 1:nind){ real pr; if (first[i] > 0) { if (last[i] > first[i]) // Recovered pr <- cell_prob(n_occasions, s[i], r[i], first[i], last[i] - 1); else // Never Recovered pr <- cell_prob(n_occasions, s[i], r[i], first[i], n_occasions); 1 ~ bernoulli(pr); } } } " ## Bundle data stan_data <- list(y = MR, nind = dim(MR)[1], n_occasions = dim(MR)[2]) ## Initial values inits <- function(){list(mean_s = runif(1, 0, 1), mean_r = runif(1, 0, 1))} ## Parameters monitored params <- c("mean_s", "mean_r") ## MCMC settings ni <- 2000 nt <- 1 nb <- 1000 nc <- 3 ## Stan model <- stan_model(model_code = mr_ss_code) mr_ss <- sampling(model, data = stan_data, init = inits, pars = params, chains = nc, thin = nt, iter = ni, warmup = nb, seed = 2, open_progress = FALSE)
8.3. The mark-recovery model fitted with the multinomial likelihood
8.3.1. Constant parameters
mr_mnl_code <- " data { int n_occasions; int marr[n_occasions, n_occasions + 1]; } parameters { real<lower=0,upper=1> mean_s; real<lower=0,upper=1> mean_r; } transformed parameters { vector[n_occasions] s; vector[n_occasions] r; simplex[n_occasions + 1] pr[n_occasions]; // Constraints for (t in 1:n_occasions) { s[t] <- mean_s; r[t] <- mean_r; } // Define the cell probabilities of the m-array // Main diagonal for (t in 1:n_occasions) { pr[t][t] <- (1 - s[t]) * r[t]; // Above main diagonal for (j in (t + 1):n_occasions) pr[t][j] <- prod(segment(s, t, j - t)) * (1 - s[j]) * r[j]; // Below main diagonal for (j in 1:(t - 1)) pr[t][j] <- 0; } // t // Last column: probability of non-recovery for (t in 1:n_occasions) pr[t][n_occasions + 1] <- 1 - sum(head(pr[t], n_occasions)); } model { // Priors mean_s ~ uniform(0, 1); // Prior for mean survival mean_r ~ uniform(0, 1); // Prior for mean recovery // Define the multinomial likelihood for (t in 1:n_occasions) marr[t] ~ multinomial(pr[t]); } " ## Bundle data stan_data <- list(marr = marr, n_occasions = dim(marr)[2] - 1) ## Initial values inits <- function(){list(mean_s = runif(1, 0, 1), mean_r = runif(1, 0, 1))} ## Parameters monitored params <- c("mean_s", "mean_r") ## MCMC settings ni <- 3000 nt <- 2 nb <- 1000 nc <- 4 ## Stan model <- stan_model(model_code = mr_mnl_code) mr <- sampling(model, data = stan_data, init = inits, pars = params, chains = nc, thin = nt, iter = ni, warmup = nb, seed = 2, open_progress = FALSE)
8.3.2. Age-dependent parameters
mr_mnl_age_code <- " data { int n_occasions; int marr_j[n_occasions, n_occasions + 1]; int marr_a[n_occasions, n_occasions + 1]; } parameters { real<lower=0,upper=1> mean_sj; real<lower=0,upper=1> mean_sa; real<lower=0,upper=1> mean_rj; real<lower=0,upper=1> mean_ra; } transformed parameters { vector[n_occasions] sj; vector[n_occasions] sa; vector[n_occasions] rj; vector[n_occasions] ra; simplex[n_occasions + 1] pr_a[n_occasions]; simplex[n_occasions + 1] pr_j[n_occasions]; // Constraints for (t in 1:n_occasions) { sj[t] <- mean_sj; sa[t] <- mean_sa; rj[t] <- mean_rj; ra[t] <- mean_ra; } // Define the cell probabilities of the juvenile m-array // Main diagonal for (t in 1:n_occasions) { pr_j[t, t] <- (1 - sj[t]) * rj[t]; // Further above main diagonal for (j in (t + 2):n_occasions) pr_j[t, j] <- sj[t] * prod(segment(sa, t + 1, j - t - 1)) * (1 - sa[j]) * ra[j]; // Below main diagonal for (j in 1:(t - 1)) pr_j[t, j] <- 0; } //t for (t in 1:(n_occasions - 1)) // One above main diagonal pr_j[t, t + 1] <- sj[t] * (1 - sa[t + 1]) * ra[t + 1]; // Last column: probability of non-recovery for (t in 1:n_occasions) pr_j[t, n_occasions + 1] <- 1 - sum(head(pr_j[t], n_occasions)); // Define the cell probabilities of the adult m-array // Main diagonal for (t in 1:n_occasions) { pr_a[t,t] <- (1 - sa[t]) * ra[t]; // Above main diagonal for (j in (t + 1):n_occasions) pr_a[t, j] <- prod(segment(sa, t, j - t)) * (1 - sa[j]) * ra[j]; // Below main diagonal for (j in 1:(t - 1)) pr_a[t, j] <- 0; } //t // Last column: probability of non-recovery for (t in 1:n_occasions) pr_a[t, n_occasions + 1] <- 1 - sum(head(pr_a[t], n_occasions)); } model { // Priors mean_sj ~ uniform(0, 1); // Prior for mean juv. survival mean_sa ~ uniform(0, 1); // Prior for mean ad. survival mean_rj ~ uniform(0, 1); // Prior for mean juv. recovery mean_ra ~ uniform(0, 1); // Prior for mean ad. recovery // Define the multinomial likelihood for (t in 1:n_occasions) { increment_log_prob(multinomial_log(marr_j[t], pr_j[t])); increment_log_prob(multinomial_log(marr_a[t], pr_a[t])); } } " ## Bundle data stan_data <- list(marr_j = marr.j, marr_a = marr.a, n_occasions = dim(marr.j)[2]-1) ## Initial values inits <- function(){list(mean_sj = runif(1, 0, 1), mean_sa = runif(1, 0, 1), mean_rj = runif(1, 0, 1), mean_ra = runif(1, 0, 1))} ## Parameters monitored params <- c("mean_sj", "mean_rj", "mean_sa", "mean_ra") # MCMC settings ni <- 2000 nt <- 1 nb <- 1000 nc <- 3 ## Stan model <- stan_model(model_code = mr_mnl_age_code) mr_age <- sampling(model, data = stan_data, init = inits, pars = params, chains = nc, thin = nt, iter = ni, warmup = nb, seed = 2, open_progress = FALSE)
8.4. Real data example: age-dependent survival in Swiss red kites
mr_mnl_age3_code <- " data { int n_age; int marr_j[n_age + 1]; int marr_a[n_age + 1]; } parameters { real<lower=0,upper=1> sjuv; real<lower=0,upper=1> ssub; real<lower=0,upper=1> sad; real<lower=0,upper=1> rjuv; real<lower=0,upper=1> rad; } transformed parameters { simplex[n_age + 1] pr_a; simplex[n_age + 1] pr_j; // Define the cell probabilities of the juvenile m-array // First element pr_j[1] <- (1 - sjuv) * rjuv; // Second element pr_j[2] <- sjuv * (1 - ssub) * rad; // Third and further elements for (t in 3:n_age) pr_j[t] <- sjuv * ssub * pow(sad, (t - 3)) * (1 - sad) * rad; // Probability of non-recovery pr_j[n_age + 1] <- 1 - sum(head(pr_j, n_age)); // Define the cell probabilities of the adult m-array // All elements for (t in 1:n_age) pr_a[t] <- pow(sad, (t - 1)) * (1 - sad) * rad; // Probability of non-recovery pr_a[n_age + 1] <- 1 - sum(head(pr_a, n_age)); } model { // Priors sjuv ~ beta(4.2, 2.8); // Informative prior for juv. survival: Analysis A //sjuv ~ uniform(0, 1); // Non-informative for juv. survival prior: Analysis B ssub ~ uniform(0, 1); // Prior for subad. survival sad ~ uniform(0, 1); // Prior for ad. survival rjuv ~ uniform(0, 1); // Prior for juv. recovery rad ~ uniform(0, 1); // Prior for ad. recovery // Define the multinomial likelihood increment_log_prob(multinomial_log(marr_j, pr_j)); increment_log_prob(multinomial_log(marr_a, pr_a)); } " ## Bundle data stan_data <- list(marr_j = marray.juv, marr_a = marray.ad, n_age = length(marray.juv) - 1) ## Initial values inits <- function(){list(sjuv = runif(1, 0, 1), ssub = runif(1, 0, 1), sad = runif(1, 0, 1), rjuv = runif(1, 0, 1), rad = runif(1, 0, 1))} ## Parameters monitored params <- c("sjuv", "ssub", "sad", "rjuv", "rad") ## MCMC settings ni <- 4000 nt <- 2 nb <- 2000 nc <- 4 ## Stan model <- stan_model(model_code = mr_mnl_age3_code) rk_ageA <- sampling(model, data = stan_data, init = inits, pars = params, chains = nc, thin = nt, iter = ni, warmup = nb, seed = 1, open_progress = FALSE)
タグ:STAn
2015-12-07 21:04
nice!(1)
コメント(0)
トラックバック(1)
コメント 0