Hi folks! I’m new to Stan and this forum. I’ve recently made 4 models modified from this: https://github.com/stan-dev/example-models/blob/master/misc/ecology/mark-recapture/cjs-K.stan
I would like to evaluate these models via WAIC or LOO with the R loo
package. For the life of me I cannot figure out how to define the log likelihood within the “generated quantities” block. Can anyone help with an example that’s applicable to the cjs-K.stan file? Intuitively I imagine incorporating bernoulli_lpmf
within a for loop
, but have no idea how to execute this. Still learning obviously! Any help is appreciated.
Hi @B-Davis! Welcome to the Stan forums - always happy to see capture-recapture work here.
To use loo
, you’ll need an object log_lik
in your Stan program as described here: https://cran.r-project.org/web/packages/loo/vignettes/loo2-with-rstan.html
In this case, you can create a vector log_lik
in the transformed parameters block that contains the log probabilities of each individual capture history. Then, increment the log probability using this vector in the model block:
/**
* Cormack-Jolly-Seber Model
*
* following section 1.2.1 of:
* http://www.maths.otago.ac.nz/home/resources/theses/PhD_Matthew_Schofield.pdf
*
*/
data {
int<lower=2> K; // capture events
int<lower=0> I; // number of individuals
int<lower=0,upper=1> X[I,K]; // X[i,k]: individual i captured at k
}
transformed data {
int<lower=0,upper=K+1> first[I]; // first[i]: ind i first capture
int<lower=0,upper=K+1> last[I]; // last[i]: ind i last capture
int<lower=0,upper=I> n_captured[K]; // n_capt[k]: num aptured at k
first = rep_array(K+1,I);
last = rep_array(0,I);
for (i in 1:I) {
for (k in 1:K) {
if (X[i,k] == 1) {
if (k < first[i])
first[i] = k;
if (k > last[i])
last[i] = k;
}
}
}
n_captured = rep_array(0,K);
for (i in 1:I)
for (k in 1:K)
n_captured[k] = n_captured[k] + X[i,k];
}
parameters {
vector<lower=0,upper=1>[K-1] phi; // phi[k]: Pr[alive at k + 1 | alive at k]
vector<lower=0,upper=1>[K] p; // p[k]: Pr[capture at k]
// note: p[1] not used in model and hence not identified
}
transformed parameters {
vector<lower=0,upper=1>[K] chi; // chi[k]: Pr[no capture > k | alive at k]
vector[I] log_lik;
{
int k;
chi[K] = 1.0;
k = K - 1;
while (k > 0) {
chi[k] = (1 - phi[k]) + phi[k] * (1 - p[k+1]) * chi[k+1];
k = k - 1;
}
}
for (i in 1:I) {
log_lik[i] = 0;
if (last[i] > 0) {
for (k in (first[i]+1):last[i]) {
log_lik[i] +=log(phi[k-1]); // i survived from k-1 to k
if (X[i,k] == 1)
log_lik[i] +=log(p[k]); // i captured at k
else
log_lik[i] +=log1m(p[k]); // i not captured at k
}
log_lik[i] +=log(chi[last[i]]); // i not seen after last[i]
}
}
}
model {
target += sum(log_lik);
}
generated quantities {
// phi[K-1] and p(K) not identified, but product is
real beta;
vector<lower=0>[K] pop_hat; // population
beta = phi[K-1] * p[K];
for (k in 1:K)
pop_hat[k] = n_captured[k] / p[k];
}
1 Like
Wow. Fantastic @mbjoseph!! Thank you!
1 Like