# CJS log likelihood

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