Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
psi <- .2
theta <- .3
p <- .5
nsite <- 100
ntime <- 4
nrep <- 3
z <- rbinom(nsite, 1, psi)
alpha <- matrix(rbinom(nsite * ntime, 1, theta), nrow = nsite)
y <- array(dim = c(nsite, ntime, nrep))
for (i in 1:nsite) {
for (j in 1:ntime) {
y[i, j, ] <- rbinom(nrep, 1, z[i] * alpha[i, j] * p)
}
}
# potential combinations of alpha that can lead to all-zero capture history
alpha_potential <- expand.grid(rep(list(c(0, 1)), ntime))
stan_d <- list(nsite = nsite,
ntime = ntime,
nrep = nrep,
y = y,
n_possible = 2^ntime,
alpha_potential = alpha_potential)
m_init <- stan_model("multiscale-occupancy.stan")
m_fit <- sampling(m_init, data = stan_d)
pairs(m_fit, pars = c("psi", "theta", "p"))
data {
int<lower = 1> nsite;
int<lower = 1> ntime;
int<lower = 1> nrep;
int<lower = 0, upper = 1> y[nsite, ntime, nrep];
int<lower = 1> n_possible;
matrix<lower = 0, upper = 1>[n_possible, ntime] alpha_potential;
}
transformed data {
int<lower = 0, upper = 1> known_present[nsite];
int<lower = 0, upper = 1> known_available[nsite, ntime];
for (i in 1:nsite) {
known_present[i] = 0;
for (j in 1:ntime) {
known_available[i, j] = 0;
for (k in 1:nrep) {
if (y[i, j, k] == 1) {
known_present[i] = 1;
known_available[i, j] = 1;
}
}
}
}
}
parameters {
real<lower = 0, upper = 1> psi;
real<lower = 0, upper = 1> theta;
real<lower = 0, upper = 1> p;
}
transformed parameters {
vector[nsite] log_lik;
{
vector[ntime] tmp_lp;
matrix[n_possible, ntime] tmp_poss;
vector[n_possible + 1] sum_poss;
for (i in 1:nsite) {
if (known_present[i]) {
for (j in 1:ntime) {
if (known_available[i, j]) {
// present and available
tmp_lp[j] = log(theta) + bernoulli_lpmf(y[i, j, ] | p);
} else {
// present, possibly unavailable
tmp_lp[j] = log_sum_exp(
log(theta) + bernoulli_lpmf(y[i, j, ] | p),
log1m(theta)
);
}
}
log_lik[i] = log(psi) + sum(tmp_lp);
} else {
// could be present or absent (was never detected)
// and there are 2^ntime possible combinations
// of alpha_{i, t} that are relevant if z_i = 1
for (jj in 1:n_possible) {
for (j in 1:ntime) {
if (alpha_potential[jj, j] == 0) {
// not available
tmp_poss[jj, j] = log1m(theta);
} else {
// available but not detected
tmp_poss[jj, j] = log(theta) + bernoulli_lpmf(y[i, j, ] | p);
}
}
sum_poss[jj] = log(psi) + sum(tmp_poss[jj, ]);
}
sum_poss[n_possible + 1] = log1m(psi);
log_lik[i] = log_sum_exp(sum_poss);
}
}
}
}
model {
target += sum(log_lik);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment