functions { int first_capture(array[] int x_i) { for (k in 1:size(x_i)) if (x_i[k] > 0) return k; return 0; } int next_capture(array[] int x_i, int now, int last) { for (k in (now + 1):last) if (x_i[k] > 0) return k; return 0; } int last_capture(array[] int x_i) { for (k_rev in 0:(size(x_i) - 1)) { int k; k = size(x_i) - k_rev; if (x_i[k] > 0) return k; } return 0; } row_vector arrival_probs_ln(int T, real mu, real sigma){ vector[T ] base; vector[T] t = linspaced_vector(T, 0.5, T - 0.5); base = -square(log(t) - mu) / (2 * square(sigma)) - log(t); return(softmax(base)'); } } data{ int N; int G; int T; array[G] int rls_day; array[N] int rls_grp; array[N, 4] int ch; int P_phi_0; int P_mu_0; int P_p_1; int P_phi_1; int P_mu_1; int P_p_2; int P_phi_2; int P_mu_2; int P_p_3; int P_phi_3; int P_mu_3; int P_p_4; // array of row vectors because matrices are stored in column-major order // so this should be faster to access row vectors in memory array[G] row_vector[P_phi_0] X_phi_0; array[G] row_vector[P_mu_0] X_mu_0; array[T] row_vector[P_p_1] X_p_1; array[T] row_vector[P_phi_1] X_phi_1; array[T] row_vector[P_mu_1] X_mu_1; array[T] row_vector[P_p_2] X_p_2; array[T] row_vector[P_phi_2] X_phi_2; array[T] row_vector[P_mu_2] X_mu_2; array[T] row_vector[P_p_3] X_p_3; array[T] row_vector[P_phi_3] X_phi_3; array[T] row_vector[P_mu_3] X_mu_3; array[T] row_vector[P_p_4] X_p_4; } transformed data { array[4, G] vector[T] M_0 = rep_array(rep_vector(0, T), 4, G); // counts of first recapture array[6, T] vector[T] M = rep_array(rep_vector(0, T), 6, T); // counts of subsequent recaptures vector[G] L_0 = rep_vector(0, G); // number never seen after release array[4] vector[T] L = rep_array(rep_vector(0, T), 4); // number never seen after capture occasion K array[4, 4] int M_i = rep_array(0, 4, 4); M_i[1, 2] = 1; //cap_hist[, 1]> 0 & cap_hist[, 2] > 0 = M[1]; M_i[2, 3] = 2; //cap_hist[, 2]> 0 & cap_hist[, 3] > 0 = M[2]; M_i[3, 4] = 3; //cap_hist[, 3]> 0 & cap_hist[, 4] > 0 = M[3]; M_i[1, 3] = 4; //cap_hist[, 1]> 0 & cap_hist[, 2] == 0 & cap_hist[, 3] > 0 = M[4]; M_i[2, 4] = 5; //cap_hist[, 2]> 0 & cap_hist[, 3] == 0 & cap_hist[, 4] > 0 = M[5]; M_i[1, 4] = 6; //cap_hist[, 1]> 0 & cap_hist[, 2:3] == 0 & cap_hist[, 4] > 0 = M[6]; for (i in 1:N){ int g; int f; int l; g = rls_grp[i]; f = first_capture(ch[i,]); l = last_capture(ch[i,]); if (f == l){ if (l == 0){ L_0[g] += 1; } else { M_0[f, g, ch[i, f]] += 1; L[ l, ch[i, l]] += 1; } } else { int now; int nxt; M_0[f, g, ch[i, f]] += 1; now = f; nxt = next_capture(ch[i,], now, l); while(nxt < l){ M[M_i[now, nxt], ch[ i, now], ch[ i, nxt]] += 1; now = nxt; nxt = next_capture(ch[i,], now, l); } M[M_i[now, l], ch[ i, now], ch[ i, l]] += 1; L[l, ch[i, l]] += 1; } } } parameters { vector[P_phi_0] Beta_phi_0; // regression parameters detection: survival 0 vector[P_phi_1] Beta_phi_1; // regression parameters detection: survival 1 vector[P_phi_2] Beta_phi_2; // regression parameters detection: survival 2 vector[P_phi_3] Beta_phi_3; // regression parameters detection: survival 3 vector[P_mu_0] Beta_mu_0; // regression parameters detection: log mean travel time vector[P_mu_1] Beta_mu_1; // regression parameters detection: log mean travel time vector[P_mu_2] Beta_mu_2; // regression parameters detection: log mean travel time vector[P_mu_3] Beta_mu_3; // regression parameters detection: log mean travel time array[4] real sigma_tt; vector[P_p_1] Beta_p_1; // regression parameters detection: collection efficiency vector[P_p_2] Beta_p_2; // regression parameters detection: collection efficiency vector[P_p_3] Beta_p_3; // regression parameters detection: collection efficiency // vector[P_p_4] Beta_p_4; // regression parameters detection: collection efficiency } transformed parameters{ vector[G] phi_0; vector[G] mu_0; // array[4] vector[T] p; array[3] vector[T] p; array[3] vector[T] phi; array[3] vector[T] mu; matrix[G, T] alpha_0 = rep_matrix(0, G, T); array[3] matrix[T, T] alpha = rep_array(rep_matrix(0, T, T), 3); for (g in 1:G){ phi_0[g] = inv_logit(X_phi_0[g] * Beta_phi_0); mu_0[g] = X_mu_0[g] * Beta_mu_0; alpha_0[g, rls_day[g]:T] = arrival_probs_ln(T - rls_day[g] + 1, mu_0[g], sigma_tt[1]); } for (t in 1:T){ phi[1, t] = inv_logit(X_phi_1[t] * Beta_phi_1); phi[2, t] = inv_logit(X_phi_2[t] * Beta_phi_2); phi[3, t] = inv_logit(X_phi_3[t] * Beta_phi_3); mu[1, t] = X_mu_1[t] * Beta_mu_1; mu[2, t] = X_mu_2[t] * Beta_mu_2; mu[3, t] = X_mu_3[t] * Beta_mu_3; p[1, t] = inv_logit(X_p_1[t] * Beta_p_1); p[2, t] = inv_logit(X_p_2[t] * Beta_p_2); p[3, t] = inv_logit(X_p_3[t] * Beta_p_3); // p[4, t] = inv_logit(X_p_4[t] * Beta_p_4); for (k in 1:3) alpha[k, t, t:T] = arrival_probs_ln(T - t + 1, mu[k, t], sigma_tt[k + 1]); } } model { array[4] matrix[G, T] Lambda_0; array[6] matrix[T, T] Lambda; array[3] vector[T] Chi; vector[G] Chi_0; //priors sigma_tt ~ normal(0, 1); Beta_phi_0[1] ~ student_t(7, 0, 2); if(P_phi_0 > 1) Beta_phi_0[2:P_phi_0] ~ student_t(7, 0, 1); Beta_mu_0[1] ~ normal(2, 1); if(P_mu_0 > 1) Beta_mu_0[2:P_mu_0] ~ normal(0, 1); Beta_p_1[1] ~ student_t(7, 0, 2); if(P_p_1 > 1) Beta_p_1[2:P_p_1] ~ student_t(7, 0, 1); Beta_phi_1[1] ~ student_t(7, 0, 2); if(P_phi_1 > 1) Beta_phi_1[2:P_phi_1] ~ student_t(7, 0, 1); Beta_mu_1[1] ~ normal(2, 1); if(P_mu_1 > 1) Beta_mu_1[2:P_mu_1] ~ normal(0, 1); Beta_p_2[1] ~ student_t(7, 0, 2); if(P_p_2 > 1) Beta_p_2[2:P_p_2] ~ student_t(7, 0, 1); Beta_phi_2[1] ~ student_t(7, 0, 2); if(P_phi_2 > 1) Beta_phi_2[2:P_phi_2] ~ student_t(7, 0, 1); Beta_mu_2[1] ~ normal(2, 1); if(P_mu_2 > 1) Beta_mu_2[2:P_mu_2] ~ normal(0, 1); Beta_p_3[1] ~ student_t(7, 0, 2); if(P_p_3 > 1) Beta_p_3[2:P_p_3] ~ student_t(7, 0, 1); Beta_phi_3[1] ~ student_t(7, 0, 2); if(P_phi_3 > 1) Beta_phi_3[2:P_phi_3] ~ student_t(7, 0, 1); Beta_mu_3[1] ~ normal(2, 1); if(P_mu_3 > 1) Beta_mu_3[2:P_mu_3] ~ normal(0, 1); // Beta_p_4[1] ~ student_t(7, 0, 2); // if(P_p_4 > 1) // Beta_p_4[2:P_p_4] ~ student_t(7, 0, 1); Chi[3] = (1 - phi[3]); Chi[2] = (1 - phi[2]) + phi[2] .* (alpha[2] * ((1 - p[3]) .* Chi[3])); Chi[1] = (1 - phi[1]) + phi[1] .* (alpha[1] * ((1 - p[2]) .* Chi[2])); Chi_0 = (1 - phi_0) + phi_0 .* (alpha_0 * ((1 - p[1]) .* Chi[1])); Lambda_0[1] = diag_pre_multiply(phi_0, alpha_0); for (k in 1:3) Lambda_0[k + 1] = diag_post_multiply(Lambda_0[k], (1 - p[k]) .* phi[k]) * alpha[k]; Lambda[1] = diag_pre_multiply(phi[1], alpha[1]); Lambda[2] = diag_pre_multiply(phi[2], alpha[2]); Lambda[3] = diag_pre_multiply(phi[3], alpha[3]); Lambda[4] = diag_post_multiply(Lambda[1], (1 - p[2]) .* phi[2]) * alpha[2]; Lambda[5] = diag_post_multiply(Lambda[2], (1 - p[3]) .* phi[3]) * alpha[3]; Lambda[6] = diag_post_multiply(Lambda[4], (1 - p[3]) .* phi[3]) * alpha[3]; target += log(Chi_0)' * L_0; for (g in 1:G){ int s = rls_day[g]; for (k in 1:3) target += (log(Lambda_0[k, g, s:]) + log(p[k, s:])') * M_0[k, g, s:]; target += log(Lambda_0[4, g, s:]) * M_0[4, g, s:]; } for (k in 1:3) target += log(Chi[k])' * L[k]; for (s in 1:T){ target += (log(Lambda[1, s, s:]) + log(p[2, s:])') * M[1, s, s:]; target += (log(Lambda[2, s, s:]) + log(p[3, s:])') * M[2, s, s:]; target += (log(Lambda[4, s, s:]) + log(p[3, s:])') * M[4, s, s:]; target += log(Lambda[3, s, s:]) * M[3, s, s:]; target += log(Lambda[5, s, s:]) * M[5, s, s:]; target += log(Lambda[6, s, s:]) * M[6, s, s:]; } }