Jolly-Seber model, state-space representation

[Peter Smits]

Hi list,

I'm trying to code the state-space representation of the JS model for capture-mark-recapture data in stan and I think I've reached the (current) limit of my mathematical ability with regards to marginalization of the (latent) discrete parameter.

The basic model is represented by three Bernoulli sampling statements.

i corresponds to observed individual, t is the point in time, p is sampling/observation probability, y_{i,t} is the observed state (0 or 1), z is the latent discrete state (0 or 1), phi is the probability of survival from t to t+1, and gamma is the probability of being "born" (recruitment).  There are no "zombies" (individuals can't be reborn).

y_{i,t} ~ Bernoulli(p_{t} * z_{i,t})

z_{i,t} ~ Bernoulli(z_{i,t-1} * phi_{t-1} + (1 - z{i,t-1}) * gamma_{t})

z_{i,1} ~ Bernoulli(gamma_{1}).

Note that I can also write z_{i,t} ~ Bernoulli(gamma_{t} + (phi_{t} - gamma{t}) * z_{i,t}) instead of the second sampling statement following Royle and Dorazio's "Hierarchical modeling and inference in ecology" book. I've also attached a relevant 2012 paper from King titled "A review of Bayesian state-space modelling of capture–recapture–recovery data". The most relevant section is 2.2 (Estimating abundance).

Given that I want to implement this in stan, I need to marginalize out the discrete parameter z. This is where me problem comes in: I don't know how to do this.

I've read over the CJS model implemented in the stan manual (page 111), and it has been quite helpful for better understanding what I'd probably have to do. I'm just stuck on how to modify that to end with the JS model via the introduction of recruitment.

Any advice would be highly appreciated.

Cheers,

Peter Smits

[Krzysztof Sakrejda]

The summation over 'z' in CJS runs over all possible unobserved outcomes after the final observation.  Similarly the summation in JS would need to run over all possible unobserved outcomes from the first possible birth up until tagging... so for one individual unobserved at t=1,2,3 and observed at t=4:

gamma(t_1) * phi(t_1) * phi(t_2) * phi(t_3) + gamma(t_2) * phi(t_2) * phi(t_3) + gamma(t_3) * phi(t_3) + gamma(t_4)

if the individual is available for recapture after birth, this becomes

gamma(t_1) * (1-p(t_1))*phi(t_1) * (1-p(t_2))* phi(t_2) * (1-p(t_3))*phi(t_3) + gamma(t_2) * (1-p(t_2))*phi(t_2) * (1-p(t_3))*phi(t_3) + gamma(t_3) * (1-p(t_3))*phi(t_3) + gamma(t_4)

That's "born at time 1 and survived 1 to 2, 2 to 3, and 3 to 4" plus "born at time 2 and survived 2 to 3 and 3 to 4" plus... etc... In both cases parts of the calculation repeat, so you can do a running calculation of the "survive and not capture" running backwards in time, keep it on the log scale, and just add the gamma component to each one, something like:

for ( t in 3:1) {
   A <- log( phi(t) ) + log(1-p(t)) + A
  B[t] <- A[t] +  log(gamma(t))
}
increment_log_prob(sum(B) + gamma(t))

It's late and I didn't check any of that but that's the gist of it.

Krzysztof


 

[Peter Smits]

Krzysztof,

Thanks for the reply. Unpacking z makes a lot of sense, thank you for pointing that out. I just want to make sure I've got this right, as I had to make an assumption regarding how to actually include my observed data.

Here is a brief model block. The full stan file is attached. 

T is total number of time windows (31 fyi), N is total number of individuals (~2000), sight is the NxT matrix of presence/absence, p is sampling, phi is survival, and gamma is birth:

model {

  real running;

  vector[T] total;

 

  running <- 0;

  for(n in 1:N) {

    for(t in T:1) {

      running <- log(phi[t]) + log(1 - p[t]) + running;

      total[t] <- running + log(gamma[t]);

      increment_log_prob(bernoulli_log(sight[n, t], 

        log_sum_exp(log(p[t]), running)));

    }

  }

}

Does this look appropriate? I'm pretty sure I messed up the increment statement, though the attached stan file does compile.

I really appreciate this, thank you.

Cheers,

Peter

[Krzysztof Sakrejda]

On Wednesday, May 27, 2015 at 11:54:29 AM UTC-4, Peter Smits wrote:

Krzysztof,

Thanks for the reply. Unpacking z makes a lot of sense, thank you for pointing that out. I just want to make sure I've got this right, as I had to make an assumption regarding how to actually include my observed data.

You're welcome, I'd love to see a good JS implementation in Stan but I have no motivation to implement one.
 

Here is a brief model block. The full stan file is attached. 

T is total number of time windows (31 fyi), N is total number of individuals (~2000), sight is the NxT matrix of presence/absence, p is sampling, phi is survival, and gamma is birth:

That's pretty big for mark-recapture.  Is this one of the big datasets that's been worked over already or a new one?

 

model {

  real running;

  vector[T] total;

 

  for(n in 1:N) {

  running <- 0;

You need to zero"running" inside this loop, otherwise this'' running'' accumulates over individuals.  If the model is not heterogeneous among individuals (only histories are) then running/total can be calculated once and put into a vector to be applied to the bernoulli_log for each individual.   Or leave it as is and make it heterogeneous later. You have enough data for that. 

    for(t in T:1) {

      running <- log(phi[t]) + log(1 - p[t]) + running;

      total[t] <- running + log(gamma[t]);

      increment_log_prob(bernoulli_log(sight[n, t],  log_sum_exp(log(p[t]), running)));

I think this can just be:
      increment_log_prob(bernoulli_log(sight[n, t],  log(p[t]) + running));

    }

  }

}

Does this look appropriate? I'm pretty sure I messed up the increment statement, though the attached stan file does compile.

That looks close, though in your place I would 1) make sure the time indexing is exactly what you want, i.e.-is phi[t] for surviving from t-1 to t or from t to t+1, and 2) Stan has a "print" statement so I would do the calculation in R for a few individuals and make sure printing the bernoulli parameter gives you what you expect. 
 

I really appreciate this, thank you.

You're welcome, if you get a reproducible example working on fake data I (and other people, I'm sure) would be interested in seeing it on the list).

Krzysztof
 

[Peter Smits]

Krzysztof,

I went back and read over Royle and Dorazio's book, and given what you've written and what I've pried out of that here is where I am at. I'm using a time forward form of the model and I've created "k" which is the related to recruitability. I've also, for the sake of simplicity, assumed time constant probabilities.

z_{i, t} = phi * z_{i, t - 1} + gamma * k

where k = prod_{k = 1}^{t - 1} (1 - z_{i, k}).

Then I've just done y_{i, t} ~ Bern(p * z_{i, t}), where p is still preservation probability.

Assuming I've done this right, I run the simulated data set through the model. What comes out the other side, however, is super unconstrained. Generally, the true values fall in the credible interval. So I think this is right? Thoughts?

Attached at a stan file, an R simulation script, and the bash script I use to run the model. Beyond the vagaries of directory trees, these should work. 

My data set is a set of global fossil occurrences from the Paleozoic. Each time point is a geologic stage and each individual is a brachiopod genus. The data set is available for free via the Paleobiology Database. Systems like brachiopods are great because they have higher preservation rate/potential than most other groups of shelly marine invertebrates. It is a non-standard data set for ecology, but fairly normal for paleobiology.

Cheers,

Peter

[Bob Carpenter]

I don't think the final model captured what you wanted
because it's not in the form of a marginalization.

More below.

> On May 26, 2015, at 4:54 PM, Peter Smits <thesti...@gmail.com> wrote:
>
> Hi list,
>
> I'm trying to code the state-space representation of the JS model for capture-mark-recapture data in stan and I think I've reached the (current) limit of my mathematical ability with regards to marginalization of the (latent) discrete parameter.
>
> The basic model is represented by three Bernoulli sampling statements.
>
> i corresponds to observed individual, t is the point in time, p is sampling/observation probability, y_{i,t} is the observed state (0 or 1), z is the latent discrete state (0 or 1), phi is the probability of survival from t to t+1, and gamma is the probability of being "born" (recruitment). There are no "zombies" (individuals can't be reborn).
>
> y_{i,t} ~ Bernoulli(p_{t} * z_{i,t})
> z_{i,t} ~ Bernoulli(z_{i,t-1} * phi_{t-1} + (1 - z{i,t-1}) * gamma_{t})
> z_{i,1} ~ Bernoulli(gamma_{1}).

Here's the way to think about this:

p(y, z, phi, p, gamma)

= p(phi, p, gamma) // prior

* PROD_{i in 1:N} bern(z[i,1] | gamma[1])

* PROD_{t in 2:T} bern(z[i,t] | z[i,t-1] * phi[t-1] + (1 - z[i, t-1]) * gamma[t])

* bern(y[i,t] | p[t] * z[i,t])

Then the goal is to marginalize out the the parameter z[i,t] representing if
individual i is alive at time t to get:

p(y, phi, p, gamma)

It's both tricky and tractable because of the constraints, which I would
think should ensure that z[i,] must be a sequence of contiguous 1s that
includes all t such that y[i,t] = 1.

* if y[i,t] = 1 then z[i,t] = 1

* if z[i,t] = 0 and any z[i,t-1] = 1, then z[i, t + n ] = 0 for all n > 0

In words, z[i,] has to be a sequence of contiguous 1s that includes
all t such that y[i,t] = 1.

So what you need to do is add all the probability for the legal z configurations
together (of which there are luckily at most T^2 due to constraints).

This is quite a bit trickier to code than the one in the manual, but has the
same principle. You want to define a function

// compute log likelihood of y for a single individual
real state_space_log(int[] y, vector phi, vector p, vector gamma) {
i.e., big sum (using log-sum-exp) over legal z values
vector[first_sighting[i] * (T - last_sighting[i])] lp;
int i;
i <- 0;
for (t_first_alive in 1:first_siting[i]) {
for (t_last_alive in (last_siting[i] + 1):T) {
log_prob[i] <- log prob for config, which involves more looping,
but now z[i,t] == (t >= t_first_alive && t <= t_last_alive)
i <- i + 1;
return log_sum_exp(log_prob);
}

It'll probably help to create an auxliary data structure of first and last
sightings, because you know it has to be all 1s before that, so then there's
just the choice of when the individual was born (which can be before the
first time point --- that's gamma[1], right?) and when the individual died (which
can be after the last time point).

And then the model is super simple.

for (i in 1:I)
y[i] ~ state_space(phi, p, gamma);

Sorry it's such a hassle! If there's only 30 time periods and 1000 individuals, it
should be tractable enough to sample.

If this isn't enough of a hint, I can help out more. I'm going to send
a different message with my implementation of the Dorazio-Royle species-site
occupancy model, which involves another marginalization. I really want
to see all these ecology models working.

- Bob

> Note that I can also write z_{i,t} ~ Bernoulli(gamma_{t} + (phi_{t} - gamma{t}) * z_{i,t}) instead of the second sampling statement following Royle and Dorazio's "Hierarchical modeling and inference in ecology" book. I've also attached a relevant 2012 paper from King titled "A review of Bayesian state-space modelling of capture–recapture–recovery data". The most relevant section is 2.2 (Estimating abundance).
>
> Given that I want to implement this in stan, I need to marginalize out the discrete parameter z. This is where me problem comes in: I don't know how to do this.
>
> I've read over the CJS model implemented in the stan manual (page 111), and it has been quite helpful for better understanding what I'd probably have to do. I'm just stuck on how to modify that to end with the JS model via the introduction of recruitment.
>
> Any advice would be highly appreciated.
>
> Cheers,
>
> Peter Smits
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <king_statespace.pdf>

[Peter Smits]

Hi Bob,

Thanks for the comments. After contemplating the model last night, you're total right in saying the previous model isn't capturing what I thought it was capturing.

Let me see if I'm understanding you're suggestion for how to implement this. I've included the first_capture and last_capture functions from the manual in order to get that information. The key is then estimating the probability of z being 1 during times NOT captured in the minimum observed duration, correct?

Through example: this observed capture history, 0000111001111000, corresponds to the minimum observed duration of 0000111111111000. Thus, the probabilities z = 1 of the 4 leading and trailing time points need to be calculated, correct?

This would lead to something similar the prob_uncaptured in the manual example of the CJS model but now with the gamma term.

Let's see if I get this right:

For observations before to the first observation 

for y[i, t=1] = 0: Pr(z[i, t] = 1) = k[t=1] = gamma[t=1] * p[t=1]

for y[i, t > 1] = 0: Pr(z[i, t > 1]) = k[t] = gamma[t] * p[t] + p[t] * (1 - phi[t - 1]) * k[t - 1]

Given that things can't be re-born/become zombies, for observations after the first observation it is possible to just use the prob_uncaptured function from the manual, right? I'm suspicious I've forgotten something key...

Honestly, I'd appreciate all the help I can get. This list and its archives have been super helpful in making Stan super accessible. I've been able to implement a lot of models, but this is the first time I've had to deal with marginalization (I don't come from a math/stats background).

Cheers,

Peter

[Bob Carpenter]

The key is this simple relation:

p(a) = SUM_b p(a,b)

To make a simpler case, let's say you have 5 observation times
and have observed the animal at time=2 and time=3, that is,

observed captures y: 01100

Then you get the following options for z

01100 first_alive=2 last_alive=3
01110
01111
11100
11110
11111 first_alive=1 last_alive=5

And what you need to do is compute p(y,z) for each of these z
on the log scale in an array,

lp[1] = log p(y=01100, z=01100 | p, phi, gamma)
lp[2] = log p(y=01100, z=01110 | p, phi, gamma)
...
lp[6] = log p(y=01100, z=11111 | p, phi, gamma)

then use:

increment_log_prob(log_sum_exp(lp));

because

log_sum_exp(lp) = log (exp(lp[1]) + exp(lp[2]) + ... + exp(lp[6]))

= log (p[1] + ... p[6])

where p[1] = exp(lp[1]) is the probability of y and the first config of z.

So you take the sum of the probabilities of each pair of y and legal z
and add them up then take the log. It's just log_sum_exp does it more
efficiently and has more arithmetic stability and accuracy.

The trick is that you're going to need to do this in a double loop over
first_alive and last_alive.

I'm busy right now trying to get the latest RStan patched for Windows,
but when that's done, I'll implement the model and do some data
simulation for this. It's not that difficult algebraically, but it's a lot of fiddly
coding with indexes.

- Bob

[Peter Smits]

Hi Bob,

I think I'm starting to get this. Would the state_space(...) function then look like (attached also in a stan file)? I apologize for poor code.

  real state_space_log(int[] y, vector phi, vector p, vector gamma) {

    vector[first_capture(y) * (size(y) - last_capture(y))] lp;

    int i;

    i <- 1;

    

    for(t_first_alive in 1:first_capture(y)) {

      for (t_last_alive in (last_capture(y) + 1):size(y)) {

        vector[2 * size(y)] sl;

        int z[size(y)];

        for(l in 1:size(y)) {

          z[l] <- 0;

        }

        for(a in t_first_alive:t_last_alive) 

          z[a] <- 1;

        

        sl[1] <- bernoulli_log(z[1], gamma[1]);

        for(j in 2:size(y)) {

          sl[j] <- bernoulli_log(z[j], z[j - 1] * (1 - phi[j - 1]) + 

                                       gamma[j] * (1 - z[j - 1]));

        }

        for(k in (size(y) + 1):(2 * size(y))) {

          sl[k] <- bernoulli_log(y[k], z[k] * p[k]);

        }

        lp[i] <- log_sum_exp(sl);

        i <- i + 1;

      }

    }

    return log_sum_exp(lp);

  }

Cheers,

Peter

[Bob Carpenter]

I think that's getting closer, but I still don't think it's
quite right. The z and y should be added together in the
same component and there should only be one log_sum_exp, where
the elements are the log likelihoods for each possible configuration
of hidden state of aliveness that is consistent with captures.

You also want to assign and save first_capture and last_capture
rather than recomputing.

I'll try to work out the full example tomorrow.

I took a look at that Ruth King paper and it's great. I love
it when domain scientists write clearly enough from a basic enough
level on the science side that I can understand what's going on.
She talks about the marginalization required to generate the likelihood
function, which is what you'd typically need for optimization (for [penalized]
MLE estimation).

- Bob

> <turnover_simple.stan>

[Bob Carpenter]

The answer I gave you before was right, but inefficient.

Reading Ruth King's paper made me think about the whole model
as an HMM and thus realize that the value chi[t] computed in
the Stan Manual's CJS model for the probability of being uncaptured (not sighted)
after time t was just the backward probability from the forward-backward
algorithm. By analogy, I'm going to define a forward variable alpha[t]
that represents the probability of being unsighted up to time t-1 but being
alive at t.

Then, given first[i] as first sighting of i and last[i] as last sighting
of i, we can use the forward value up to first[i] and backward value back
to last[i] and condition everything in the middle knowing that all the x[i,t] = 1
(the individual can't die and come back to life between first[i] and last[i]).

So I'd probably write this as a likelihood function for y conditioned on
a whole bunch of things:

real
log_likelihood(int[] y, // observations
int first, // first sighting time
int last, // last sighting time
vector alpha, // in study at first sighting
vector chi, // unsighted after last sighting
vector phi, // survive from t to t+1
vector p) { // sighted given alive
real ll;
ll <- alpha[first] + chi[last];
for (t in first:(last-1))
ll <- ll + bernoulli_log(1, phi[t]);
for (t in first:last)
ll <- ll + bernoulli_log(y[t], p[t]);
return ll;
}

The first loop is that it survives from each point to the next
point between first and last (you don't know if you survive beyond
last, which is why phi[last] isn't included). The second loop just
takes the sighting/capture probabilities for each point.

The first, last, and alpha values all get computed as in the previous model.
I'd define first and last as transformed data so as not to have to recompute.

Im assuming alpha and chi are defined on the log scale.

You can then define alpha and chi as locals in the model block, where
b[t] is chance of entering study at time t.

alpha[1] <- log(b[1]);
for (t in 2:T)
alpha[t] <- log_sum_exp(alpha[t-1] + bernoulli_log(1, phi[t-1]),
b[t]);

The log_sum_exp is for the mixture of either being in the study at time t-1
and surviving to time to or entering the study at time t.

Moving on to the abundance study, you'll also have a bunch of individuals
who were never sighted. So you have to increment the log density by the
probability of that, too. Which is another calculation like the others,
because they can come into the study or exit the study, but at each point they
are in the study, they are unobserved.

I think that's it. Let me know if you have more questions.

This is fun! Especially when I don't have to debug it all. You
really want to check all the boundary conditions in the above to make
sure they're right.

- Bob

[Peter Smits]

Bob,

Thanks for the extensive help with this. I really appreciate it.

I've implemented the earlier, less efficient, suggestion. I made the changes you suggested, and I think it works/makes sense now. I'm sticking with the less efficient version for now as I currently understand it better, but I will most likely revisit efficiency later after I've done some more testing etc.

I'm glad the King paper has been useful. Of the paper's I've found on this topic, it was one of the better recent reviews.

Cheers,

Peter

[Peter Smits]

Hi again,

Attached is the stan model file I've been using. Whenever I try to use it, the model fails during initialization. The three parameters of interest are constrained to be between 0 and 1, but I'm getting a message that "Initialization between (-2, 2) failed after 100 attempts." Is there something still wrong with the sampling statement I've written?

I've also attached the shell script I use to run the (cmd)stan model and a small dataset (n = 100, t = 31).

Sorry for the continued problems.

Cheers,

Peter

[Bob Carpenter]

King also wrote the chapter on ecology models for the MCMC
handbook:

http://www.mcmchandbook.net/HandbookTableofContents.html

But I don't see a version online.

I'm going to see if I can get that ring recovery data and fit some of
the recovery models and temperature covariate models.

I'll probably get around to writing the efficient form of the
JS model in a few weeks as an example.

- Bob

[Bob Carpenter]

I'd guess one of the bernoulli_log is getting a parameter
that's incompatible with its outcome.

If you do a print() inside the function, you should be
able to track down where it's happening.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <turnover_simple.stan><turnover_model_wrap.sh><sight_info.data.R>

[Peter Smits]

Hi Bob,

Thanks for the tip. I've used print to check the values for the different variables and steps in the sampling statement. What I've learned is that there is never an illegal parameter state (always between 0 and 1), which is good. I've changed the init argument for the cmdstan call to init=1 just to help limit possible concerns. I've not really found anything obvious which could be causing this problem, though I did fix a tiny error in the math. This didn't change anything, however.

What i have noticed is that the parameter estimate aren't updating between iterations. I have no idea what could be causing this (flat likelihood surface?), but I'm guessing this is the big problem. My code will print out the vectors with for sample, and they all appear to be valid values for the parameters. But, after it has all run, I still get the initialization failed message and the model stops.

Would this be a situation where putting non-uniform priors on, say, logit(gamma) might help with inference? Any help would be appreciated.

Sample from print output.

phi = [0.360966,0.35698,0.325728,0.269452,0.581838,0.474542,0.617086,0.456718,0.627167,0.520624,0.650452,0.676364,0.679755,0.400965,0.636721,0.343335,0.345352,0.511662,0.286681,0.673628,0.57874,0.388223,0.531224,0.727612,0.575721,0.566106,0.546064,0.434811,0.417624,0.354157]

gamma = [0.618148,0.571221,0.723834,0.394825,0.272438,0.660337,0.290158,0.672112,0.582204,0.312558,0.609324,0.399239,0.606841,0.282272,0.587201,0.618179,0.388153,0.422873,0.343704,0.36921,0.535746,0.514768,0.406327,0.32313,0.675828,0.472903,0.597111,0.549852,0.694149,0.594616,0.595549]

p = [0.396419,0.550119,0.626042,0.727947,0.679115,0.50664,0.334636,0.394073,0.653771,0.304764,0.418969,0.659981,0.373409,0.390601,0.492196,0.634144,0.678902,0.31745,0.529694,0.509664,0.3234,0.533394,0.477587,0.402319,0.4981,0.667508,0.393094,0.587889,0.320695,0.490758,0.593029]

phi = [0.360966,0.35698,0.325728,0.269452,0.581838,0.474542,0.617086,0.456718,0.627167,0.520624,0.650452,0.676364,0.679755,0.400965,0.636721,0.343335,0.345352,0.511662,0.286681,0.673628,0.57874,0.388223,0.531224,0.727612,0.575721,0.566106,0.546064,0.434811,0.417624,0.354157]

gamma = [0.618148,0.571221,0.723834,0.394825,0.272438,0.660337,0.290158,0.672112,0.582204,0.312558,0.609324,0.399239,0.606841,0.282272,0.587201,0.618179,0.388153,0.422873,0.343704,0.36921,0.535746,0.514768,0.406327,0.32313,0.675828,0.472903,0.597111,0.549852,0.694149,0.594616,0.595549]

p = [0.396419,0.550119,0.626042,0.727947,0.679115,0.50664,0.334636,0.394073,0.653771,0.304764,0.418969,0.659981,0.373409,0.390601,0.492196,0.634144,0.678902,0.31745,0.529694,0.509664,0.3234,0.533394,0.477587,0.402319,0.4981,0.667508,0.393094,0.587889,0.320695,0.490758,0.593029]

phi = [0.360966,0.35698,0.325728,0.269452,0.581838,0.474542,0.617086,0.456718,0.627167,0.520624,0.650452,0.676364,0.679755,0.400965,0.636721,0.343335,0.345352,0.511662,0.286681,0.673628,0.57874,0.388223,0.531224,0.727612,0.575721,0.566106,0.546064,0.434811,0.417624,0.354157]

gamma = [0.618148,0.571221,0.723834,0.394825,0.272438,0.660337,0.290158,0.672112,0.582204,0.312558,0.609324,0.399239,0.606841,0.282272,0.587201,0.618179,0.388153,0.422873,0.343704,0.36921,0.535746,0.514768,0.406327,0.32313,0.675828,0.472903,0.597111,0.549852,0.694149,0.594616,0.595549]

p = [0.396419,0.550119,0.626042,0.727947,0.679115,0.50664,0.334636,0.394073,0.653771,0.304764,0.418969,0.659981,0.373409,0.390601,0.492196,0.634144,0.678902,0.31745,0.529694,0.509664,0.3234,0.533394,0.477587,0.402319,0.4981,0.667508,0.393094,0.587889,0.320695,0.490758,0.593029]

phi = [0.360966,0.35698,0.325728,0.269452,0.581838,0.474542,0.617086,0.456718,0.627167,0.520624,0.650452,0.676364,0.679755,0.400965,0.636721,0.343335,0.345352,0.511662,0.286681,0.673628,0.57874,0.388223,0.531224,0.727612,0.575721,0.566106,0.546064,0.434811,0.417624,0.354157]

gamma = [0.618148,0.571221,0.723834,0.394825,0.272438,0.660337,0.290158,0.672112,0.582204,0.312558,0.609324,0.399239,0.606841,0.282272,0.587201,0.618179,0.388153,0.422873,0.343704,0.36921,0.535746,0.514768,0.406327,0.32313,0.675828,0.472903,0.597111,0.549852,0.694149,0.594616,0.595549]

p = [0.396419,0.550119,0.626042,0.727947,0.679115,0.50664,0.334636,0.394073,0.653771,0.304764,0.418969,0.659981,0.373409,0.390601,0.492196,0.634144,0.678902,0.31745,0.529694,0.509664,0.3234,0.533394,0.477587,0.402319,0.4981,0.667508,0.393094,0.587889,0.320695,0.490758,0.593029]

phi = [0.360966,0.35698,0.325728,0.269452,0.581838,0.474542,0.617086,0.456718,0.627167,0.520624,0.650452,0.676364,0.679755,0.400965,0.636721,0.343335,0.345352,0.511662,0.286681,0.673628,0.57874,0.388223,0.531224,0.727612,0.575721,0.566106,0.546064,0.434811,0.417624,0.354157]

gamma = [0.618148,0.571221,0.723834,0.394825,0.272438,0.660337,0.290158,0.672112,0.582204,0.312558,0.609324,0.399239,0.606841,0.282272,0.587201,0.618179,0.388153,0.422873,0.343704,0.36921,0.535746,0.514768,0.406327,0.32313,0.675828,0.472903,0.597111,0.549852,0.694149,0.594616,0.595549]

p =[0.396419,0.550119,0.626042,0.727947,0.679115,0.50664,0.334636,0.394073,0.653771,0.304764,0.418969,0.659981,0.373409,0.390601,0.492196,0.634144,0.678902,0.31745,0.529694,0.509664,0.3234,0.533394,0.477587,0.402319,0.4981,0.667508,0.393094,0.587889,0.320695,0.490758,0.593029]

[Bob Carpenter]

The probably cause is a rejection due to numerical issues
somewhere. Given all the work's being done in one function,
that's likely the source of problem.

Look at the return value and make sure it's finite. If
it is, then I don't see what the problem could be.

- Bob

[Peter Smits]

You're right. Looking log_sum_exp(lp) where lp is a vector of log probabilities, it sometimes has a value of -inf. 

I think the issue is the definition of lp. it is defined vector[first_capture(y) * (size(y) - last_capture(y))] lp; What happens, is that if last_capture is equal to size(y) the vector has length 0. And that is bad. This causes further problems when trying to loop from (last_capture + 1):size(y). 

So, I add 1 to the length of the vector lp, then it is no longer degenerate in the case where last_capture(y) == size(y). I then also changef the t_last_alive loop to last_capture:S. This seems to have fixed it.

The stan file is attached.

Thanks again, Bob. Your help has been really illuminating.

Cheers,

Peter

[Bob Carpenter]

As long as there's no undefined variables floating around.

I'd think what you want to do is this:

vector[first_capture(y) * max(size(y) - last_capture(y), 1)] ...;

BUT, that runs into a bug we have in using max() within bounds
that I'm currently in the process of fixing. So you'll have to do it like
this:

vector[foo(first_capture(y), size(y), last_capture(y))] ...;

with

int foo(int fc, int s, int lc) {
return fc * max(1, s - lc);
}

The fix for max in vector declarations probably won't make it into
Stan 2.7.0, which we're trying to rush out the door to deal
with BH (Boost headers) changing versions on us so that nobody can
get RStan to build.

- Bob