Hey again,
I was able to compute the likelihood for each individual without zeros. I did this by only computing the marginal likelihood of the alive states until the last primary that individuals were observed. When individuals were not observed in a secondary, the diagnostic (data) process was excluded and only the “not captured” probability (1 - p) was propagated. Only after individuals weren’t captured was the dead state updated. I was also able to use matrix multiplication for lots of it. Here’s the code which works great in a partial sum:
real partial_sum_lpmf(data array[] int seq_ind, int start, int end,
data array[,,,] int y, data array[,,] int q, data int n_s, data int n_o, data int n_d, data int n_prim, data int n_sec, data array[] int first, data array[] int last, data array[] int first_sec, data vector tau,
real eta_a, vector phi_a, real phi_b, vector psi_a, vector p_a, vector r, vector fr, matrix m, array[] matrix n) {
// number of individuals in partial sum and initialise partial target
int n_ind = end - start + 1;
real ptarget = 0;
// parameters
real eta;
array[n_s - 1] vector[n_prim - 1] phi, psi;
array[n_s - 1] matrix[n_sec, n_prim] p;
tuple(vector[n_prim], matrix[n_sec, n_prim]) delta;
// transform intercepts
vector[n_s - 1] log_phi_a = log(phi_a);
// transition rate/probability matrices (trm/tpm)
tuple(array[n_prim - 1] matrix[n_s, n_s],
array[n_prim, n_sec] matrix[n_s - 1, n_o],
array[n_prim, n_sec] matrix[n_o - 1, n_d - 1]) tpm;
array[n_prim - 1] matrix[n_s, n_s] trm;
// initialise probabilities
array[n_prim] matrix[n_o - 1, n_sec] acc;
matrix[n_s, n_prim] gamma;
// for each individual
for (i in 1:n_ind) {
// index that maps to original n_ind
int ii = i + start - 1;
// initialise marginal probabilities
gamma = rep_matrix(1, n_sec, n_prim);
// probability of entering as infected
eta = eta_a;
// marginal probability of alive states at first capture
gamma[1:2, first[ii]] = [ 1 - eta, eta ]';
// mortality rates
phi[1] = rep_vector(phi_a[1], n_prim - 1);
phi[2] = exp(log_phi_a[2] + phi_b * m[1:(n_prim - 1), ii]);
// infection state transition rates
psi[1] = rep_vector(psi_a[1], n_prim - 1);
psi[2] = rep_vector(psi_a[2], n_prim - 1);
// for successive primary occasions
for (t in first[ii]:(n_prim - 1)) {
// ecological trm
trm[t][1, 1] = -(psi[1][t] + phi[1][t]);
trm[t][1, 2] = psi[2][t];
trm[t][2, 1] = psi[1][t];
trm[t][2, 2] = -(psi[2][t] + phi[2][t]);
trm[t][3, 1] = phi[1][t];
trm[t][3, 2] = phi[2][t];
trm[t][:, 3] = zeros_vector(3);
// ecological tpm
tpm.1[t] = matrix_exp(trm[t] * tau[t]);
} // t
// individual and sample infection detection probabilities
delta.1 = 1 - pow(1 - r[1], m[:, ii]);
delta.2 = 1 - pow(1 - r[2], n[ii]);
// detection probabilities (fixed at 1 for secondary of first capture)
for (s in 1:(n_s - 1)) {
p[s] = rep_matrix(p_a[s], n_sec, n_prim);
p[s][first_sec[ii], first[ii]] = 1;
} // s
// for each primary occasion
for (t in first[ii]:n_prim) {
// for each secondary occasion
for (k in 1:n_sec) {
// observation tpm
tpm.2[t, k][1, 1] = p[1][k, t] * (1 - fr[1]);
tpm.2[t, k][1, 2] = p[1][k, t] * fr[1];
tpm.2[t, k][1, 3] = 1 - p[1][k, t];
tpm.2[t, k][2, 1] = p[2][k, t] * (1 - delta.1[t]);
tpm.2[t, k][2, 2] = p[2][k, t] * delta.1[t];
tpm.2[t, k][2, 3] = 1 - p[2][k, t];
// if individual was not detected
if (q[ii, t, k] == 0) {
// marginal detection probabilities
gamma[1:2, t] .*= tpm.2[t, k][:, 3];
} else {
// diagnostic tpm
tpm.3[t, k][1, 1] = 1 - fr[2];
tpm.3[t, k][1, 2] = fr[2];
tpm.3[t, k][2, 1] = 1 - delta.2[k, t];
tpm.3[t, k][2, 2] = delta.2[k, t];
// probability of each alive observed state conditioned on data
for (o in 1:(n_o - 1)) {
acc[t][o, k] = prod(tpm.3[t, k][o, y[ii, t, k]]);
} // o
// marginal probability after each secondary for each alive ecological state
gamma[1:2, t] .*= tpm.2[t, k][1:2, 1:2] * acc[t][1:2, k];
}
} // k
} // t
// marginal probability of alive states between first and last primary with captures
if (first[ii] < last[ii]) {
for (t in (first[ii] + 1):last[ii]) {
gamma[1:2, t] .*= tpm.1[t - 1][1:2, 1:2] * gamma[1:2, t - 1];
} // t
}
// if last captured in the last primary
if (last[ii] == n_prim) {
// increment target with marginal probabilities of alive states
ptarget += log(sum(gamma[1:2, n_prim]));
} else {
// marginal probabilities of all ecological states in primary after last capture
gamma[1:2, last[ii] + 1] .*= tpm.1[last[ii]][1:2, 1:2] * gamma[1:2, last[ii]];
gamma[3, last[ii] + 1] = tpm.1[last[ii]][3, 1:2] * gamma[1:2, last[ii]];
// if first occasion after last capture is the last primary
if ((last[ii] + 1) == n_prim) {
// increment target with marginal probabilities of all ecological states
ptarget += log(sum(gamma[:, n_prim]));
} else {
// marginal probabilities of all ecological states until last primary
for (t in (last[ii] + 2):n_prim) {
gamma[:, t] .*= tpm.1[t - 1] * gamma[:, t - 1];
} // t
// increment target with marginal probabilities of all ecological states
ptarget += log(sum(gamma[:, n_prim]));
}
}
} // i
return(ptarget);
}
I once again tried to translate this to log scale. As far as I can tell, I just need to sum the log probs instead of multiply, log_sum_exp() the log probs instead of sum, and write a function that does matrix multiplication for log matrices (that is, compute log(AB) from log(A) and log(B)), which is here (note I wrote two other ones that take (row_)vectors:
/**
* Compute log(A * B) from log(A) and log(B)
*
* @param a: Logarithm of matrix A
* @param b: Logarithm of matrix B
*
* @return Logarithm of matrix A * B
*/
matrix log_mat_prod(matrix a, matrix b) {
int x = rows(a);
int y = cols(b);
matrix[x, y] c;
for (i in 1:x) {
for (j in 1:y) {
c[i, j] = log_sum_exp(row(a, i)' + col(b, j));
}
}
return(c);
}
Now, it’s easy to translate the first code block to work with logs. Note that I do initialise gamma with 0s (and 1s in the first version) to make easier to iterate gamma in the observation process (i.e., gamma[1:2, t] += log(tpm.2[t, k][:, 3]); however, none of these 0 values ever get passed to (p )target, but perhaps it’s still an issue?
real partial_sum_lpmf(data array[] int seq_ind, int start, int end,
data array[,,,] int y, data array[,,] int q, data int n_s, data int n_o, data int n_d, data int n_prim, data int n_sec, data array[] int first, data array[] int last, data array[] int first_sec, data vector tau,
real eta_a, vector phi_a, real phi_b, vector psi_a, vector p_a, vector r, vector fr, matrix m, array[] matrix n) {
// number of individuals in partial sum and initialise partial target
int n_ind = end - start + 1;
real ptarget = 0;
// parameters
real eta;
array[n_s - 1] vector[n_prim - 1] phi, psi;
array[n_s - 1] matrix[n_sec, n_prim] p;
tuple(vector[n_prim], matrix[n_sec, n_prim]) delta;
// transform intercepts
vector[n_s - 1] log_phi_a = log(phi_a);
// transition rate/probability matrices (trm/tpm)
tuple(array[n_prim - 1] matrix[n_s, n_s],
array[n_prim, n_sec] matrix[n_s - 1, n_o],
array[n_prim, n_sec] matrix[n_o - 1, n_d - 1]) tpm;
array[n_prim - 1] matrix[n_s, n_s] trm;
// initialise probabilities
array[n_prim] matrix[n_o - 1, n_sec] acc;
matrix[n_s, n_prim] gamma;
// for each individual
for (i in 1:n_ind) {
// index that maps to original n_ind
int ii = i + start - 1;
// initialise marginal log probabilities
gamma = rep_matrix(0, n_sec, n_prim);
// probability of entering as infected
eta = eta_a;
// marginal log probability of alive states at first capture
gamma[1:2, first[ii]] = log([ 1 - eta, eta ]');
// mortality rates
phi[1] = rep_vector(phi_a[1], n_prim - 1);
phi[2] = exp(log_phi_a[2] + phi_b * m[1:(n_prim - 1), ii]);
// infection state transition rates
psi[1] = rep_vector(psi_a[1], n_prim - 1);
psi[2] = rep_vector(psi_a[2], n_prim - 1);
// for successive primary occasions
for (t in first[ii]:(n_prim - 1)) {
// ecological trm
trm[t][1, 1] = -(psi[1][t] + phi[1][t]);
trm[t][1, 2] = psi[2][t];
trm[t][2, 1] = psi[1][t];
trm[t][2, 2] = -(psi[2][t] + phi[2][t]);
trm[t][3, 1] = phi[1][t];
trm[t][3, 2] = phi[2][t];
trm[t][:, 3] = zeros_vector(3);
// ecological tpm
tpm.1[t] = matrix_exp(trm[t] * tau[t]);
} // t
// individual and sample infection detection probabilities
delta.1 = 1 - pow(1 - r[1], m[:, ii]);
delta.2 = 1 - pow(1 - r[2], n[ii]);
// detection probabilities (fixed at 1 for secondary of first capture)
for (s in 1:(n_s - 1)) {
p[s] = rep_matrix(p_a[s], n_sec, n_prim);
p[s][first_sec[ii], first[ii]] = 1;
} // s
// for each primary occasion
for (t in first[ii]:n_prim) {
// for each secondary occasion
for (k in 1:n_sec) {
// observation tpm
tpm.2[t, k][1, 1] = p[1][k, t] * (1 - fr[1]);
tpm.2[t, k][1, 2] = p[1][k, t] * fr[1];
tpm.2[t, k][1, 3] = 1 - p[1][k, t];
tpm.2[t, k][2, 1] = p[2][k, t] * (1 - delta.1[t]);
tpm.2[t, k][2, 2] = p[2][k, t] * delta.1[t];
tpm.2[t, k][2, 3] = 1 - p[2][k, t];
gamma[1:2, t] = zeros_vector(2);
// if individual was not detected
if (q[ii, t, k] == 0) {
// marginal detection probabilities
gamma[1:2, t] += log(tpm.2[t, k][:, 3]);
} else {
// diagnostic tpm
tpm.3[t, k][1, 1] = 1 - fr[2];
tpm.3[t, k][1, 2] = fr[2];
tpm.3[t, k][2, 1] = 1 - delta.2[k, t];
tpm.3[t, k][2, 2] = delta.2[k, t];
// probability of each alive observed state conditioned on data
for (o in 1:(n_o - 1)) {
acc[t][o, k] = sum(log(tpm.3[t, k][o, y[ii, t, k]]));
} // o
// marginal probability after each secondary for each alive ecological state
gamma[1:2, t] += log_mat_prod(log(tpm.2[t, k][1:2, 1:2]), acc[t][1:2, k]);
}
} // k
} // t
// marginal probability of alive states between first and last primary with captures
if (first[ii] < last[ii]) {
for (t in (first[ii] + 1):last[ii]) {
gamma[1:2, t] += log_mat_prod(log(tpm.1[t - 1][1:2, 1:2]), gamma[1:2, t - 1]);
} // t
}
// if last captured in the last primary
if (last[ii] == n_prim) {
// increment target with marginal probabilities of alive states
ptarget += log_sum_exp(gamma[1:2, n_prim]);
} else {
// marginal probabilities of all ecological states in primary after last capture
gamma[1:2, last[ii] + 1] += log_mat_prod(log(tpm.1[last[ii]][1:2, 1:2]), gamma[1:2, last[ii]]);
gamma[3, last[ii] + 1] = log_mat_prod(log(tpm.1[last[ii]][3, 1:2]), gamma[1:2, last[ii]]);
// if first occasion after last capture is the last primary
if ((last[ii] + 1) == n_prim) {
// increment target with marginal probabilities of all ecological states
ptarget += log_sum_exp(gamma[:, n_prim]);
} else {
// marginal probabilities of all ecological states until last primary
for (t in (last[ii] + 2):n_prim) {
gamma[:, t] += log_mat_prod(log(tpm.1[t - 1]), gamma[:, t - 1]);
} // t
// increment target with marginal probabilities of all ecological states
ptarget += log_sum_exp(gamma[:, n_prim]);
}
}
} // i
return(ptarget);
}
I guess my issue comes down to this: when can you have a 0 in gamma and when can’t you?