Estimate a model with a long vector of latent discrete parameters in Stan

[Ben Lambert]

Hi,

I am estimating a model for animal movement where in each time step the animal is in one of two Z states:

- '0' state - foraging: the animal searches for food; moving relatively short distances per time step.

- '1' state - Migrating: the animal is moving to a far-away destination, and hence moves long distances per time step.

I attach fake data for this problem (the true state of the animal is in the 'Z' part of the list, and the distance moved is in the 'l' part of the list). Visually it is fairly obvious as to which state the animal is in at any point in time. I would like to use Stan to try to estimate the state of the animal at a given time step, given its distance moved in this step and all past steps. 

I am generating the fake data by assuming that the state is persistent; an animal remains foraging or migrating for multiple time steps. I want to estimate a model in Stan that similarly allows for a persistence in State.

The issue here is that the state is discrete, and more troubling, it takes on a value for each of the time periods (I have chosen N=200 for the fake data I attach, although the real data is around 10000 time steps.) The issue with this long data is that I can't see a way to marginalise out the latent variable. Instead at the moment I have tried to use a continuous latent variable, constrained between 0 and 1, with a switch at 0.5 that dictates the actual state of the system for each time point. The code for this model is below.

Unsurprisingly this model doesn't work. The chains for the 'Z_continuous' parameters aren't mixing at all; I suppose this is because the log-probability does not vary with this variable unless it crosses 0.5.

Does anyone have an idea as to how I might be able to estimate this model, or a similarly-behaving one in Stan?

Further information: I have estimated the model in JAGS, and it works well; the model is very good at returning the state 'Z' at any point in time. The issue however is speed; JAGS is slow, and as such will be impossible to use to estimate the full data set. If I can't get this to work in Stan then I may try 'pymc3', although have never used it before. It also occurs to me that I could do some sort of particle MC sampling here. However, any other ideas are also very welcome!

Best,

Ben

data {

  int N; ## Num obs

  real<lower=0> l[N]; ## Distance moved

  real<lower=0> aPrior0;

  real<lower=0> aPriorExtra;

}

parameters{

  real<lower=0,upper=1> Z_continuous[N];

  real<lower=0> a0;

  real<lower=0> aExtra;

  real<lower=0> sigma;

}

transformed parameters{

  real Z[N];

  for(t in 1:N){

    Z[t] = Z_continuous[t] > 0.5;

  }

}

model{

  for(t in 2:N){

    if(Z_continuous[t]<0.5){

      l[t] ~ exponential(a0+aExtra);

    } else{

      l[t] ~ exponential(a0);

    }

    Z_continuous[t] ~ normal(Z_continuous[t-1],sigma);

  }

  a0 ~ exponential(aPrior0);

  aExtra ~ exponential(aPriorExtra);

}

 

[Mike Lawrence]

Have you considered a mixture model over the two states? Something along the lines of:

data {

  int N; ## Num obs

  real<lower=0> l[N]; ## Distance moved

  real<lower=0> aPrior0;

  real<lower=0> aPriorExtra;

}

parameters{

  real state_logit[N]; //mixture probs on logit scale 

  real<lower=0> a0;

  real<lower=0> aExtra;

  real<lower=0> sigma;

}

model{

  a0 ~ exponential(aPrior0);

  aExtra ~ exponential(aPriorExtra);

  state_logit[1] ~ normal(0,sigma); // you seemed to forget this, so I'm guessing at what you might want as the prior here

  for(t in 2:N){

    state_logit[t] ~ normal(state_logit[t-1],sigma);

    target += log_sum_exp(

        log(inv_logit(state_logit[t])) + exponential_lpdf( l | a0 + aExtra )

        , log1m(inv_logit(state_logit[t])) + exponential_lpdf( l | a0 )

    ) ;

  }

}

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

Good timing. I was at ISEC watching all the movement
HMM talks and decided to code it with simulated data
myself. It's just an HMM, so it's dead easy to marginalize.

I followed the conventions on change-in-bearing
being in (-pi, pi) of the moveHMM R package (on CRAN).

data {
int<lower=0> N;
vector<lower=-pi(), upper=pi()>[N] turn;
vector<lower=0>[N] dist;
int<lower=1> K;
}
parameters {
simplex[K] theta[K];
vector<lower=-pi(), upper=pi()>[K] mu;
vector<lower=0>[K] kappa;
positive_ordered[K] alpha;
vector<lower=0>[K] sigma;
}
model {
vector[K] log_theta[K];
vector[K] lp;

for (k in 1:K) log_theta[k] = log(theta[k]);
lp = rep_vector(-log(K), K);
for (n in 1:N) {
for (k in 1:K)
lp[k]
= log_sum_exp(log_theta[k] + lp)
+ von_mises_lpdf(turn[n] | mu[k], kappa[k])
+ lognormal_lpdf(dist[n] | alpha[k], sigma[k]);
}
target += log_sum_exp(lp);

kappa ~ lognormal(0, 2);
alpha ~ lognormal(0, 2);
sigma ~ lognormal(0, 1);
}

It uses the new Stan 2.10 model notation.

lp stores a column of values in the forward-backward algorithm
on the log scale.

see: https://en.wikipedia.org/wiki/Forward–backward_algorithm

You can see that through this relation (a is the transition
matrix and b is the current log densities of ending in
the particular state given the input data):

log(a' * b) = log(sum(a .* b))

= log(sum(a[1] * b[1] + ... + a[N] * b[N]))

= log(sum(exp(log(a[1] * b[1]) + ... + exp(log(a[N] * b[N])))))

= log(sum(exp(log(a[1]) + log(b[1])) + ...)

= log_sum_exp(log(a[1]) + log(b[1]), ..., log(a[N]) + log(b[N]))

= log_sum_exp(log(a) + log(b))

At the end, I only need to log_sum_exp the final forward state to
get the final likelihood for (turn[n], dist[n]) (so it's not a bug
that the target += is inot in the loop).

I didn't think too hard about priors. It does a reasonable,
but not great, job of recovering simulated parameters. And you
can swap in Weibull or whatever makes sense for the step sizes.
I just used a von Mises on the change in bearing---our von Mises
lets you choose your own interval.

And you can easily turn the transitions into multi-logit regressions.

I ran the simulation for N = 10,000 and it worked fine and ran in
about 6 minutes on my Mac notebook (using RStan 2.10) and got
a reasonable n_eff using a single chain for 200 iterations:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1,1] 0.91 0.00 0.00 0.90 0.91 0.91 0.92 0.92 30 1.00
theta[1,2] 0.09 0.00 0.00 0.08 0.08 0.09 0.09 0.10 30 1.00
theta[2,1] 0.07 0.00 0.00 0.06 0.07 0.07 0.07 0.07 61 1.01
theta[2,2] 0.93 0.00 0.00 0.93 0.93 0.93 0.93 0.94 61 1.01
mu[1] -0.21 0.01 0.05 -0.30 -0.24 -0.21 -0.18 -0.12 72 0.99
mu[2] -0.02 0.00 0.01 -0.03 -0.02 -0.02 -0.01 0.00 86 1.07
kappa[1] 0.49 0.00 0.02 0.43 0.47 0.49 0.50 0.52 50 1.00
kappa[2] 2.96 0.00 0.05 2.88 2.92 2.96 2.99 3.05 100 0.99
alpha[1] 1.00 0.00 0.00 1.00 1.00 1.00 1.01 1.01 100 0.99
alpha[2] 2.78 0.00 0.01 2.76 2.77 2.78 2.78 2.80 100 0.99
sigma[1] 0.25 0.00 0.00 0.25 0.25 0.25 0.25 0.26 100 0.99
sigma[2] 0.78 0.00 0.01 0.77 0.78 0.79 0.79 0.81 100 1.00

vs. simulated:

alpha_dist <- c(1, 3);
sigma_dist <- c(0.25, 0.5);

mu_turn <- c(-0.2, 0);
kappa_turn <- c(0.5, 4);

You'll probably want to think harder about priors and likelihood to
avoid the sampler or optimizer getting stuck in a mode (or at a likelihood
spike at a boundary). Also, this is just using
the precision-like value kappa (bigger is lower variance).

From examples I've run, it tends to overestimate
the time it stays in a state and overestimate the deviation
on the movement sizes, but I haven't done a controlled study.
I think this may be inherent to the kind of mixture modeling
being done here, with the bias coming from the 90% remain-in-state
probabilities making each state want to soak up boundary cases
it wouldn't otherwise soak up. But I'm not sure.

I'll attach the Stan file and the data simulator and the code to run it.
I made the simulated data easily separable just to make sure the algorithm
was coded properly.

I'm also attaching a file to read from the moveHMM sample data,
which also fits OK with the Weibull that they use (I didn't use
the covariates). I wanted to go through the logic of turning
angle to make sure I did it right.

I'll write this all up as a regular case study when I get the chance.

- Bob

[Ben Lambert]

Hi Bob,

That's great - thanks very much for this. The elk study (I'm sure you know) is very much the canonical paper on this topic, so very useful to see it.

Just trying to understand it and had a few questions:

1. The model you provide estimates the transition matrix for switching between the two states of the system across all time. Is there a way to estimate the state of the system at a given point in time? Essentially this would be a continuous variable of length N that represents the probability of being in one of two states.

2. Is this model extendable to allow a switching probability that depends on covariates?

Best,

Ben

[Bob Carpenter]

> On Jul 5, 2016, at 8:46 AM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Hi Bob,
>
> That's great - thanks very much for this. The elk study (I'm sure you know) is very much the canonical paper on this topic, so very useful to see it.
>
> Just trying to understand it and had a few questions:
>
> 1. The model you provide estimates the transition matrix for switching between the two states of the system across all time. Is there a way to estimate the state of the system at a given point in time? Essentially this would be a continuous variable of length N that represents the probability of being in one of two states.

Yes, you just need the backward part of the forward-backward
algorithm. I can put that in the case study. It's just as
straightforward as the forward algorithm.

> 2. Is this model extendable to allow a switching probability that depends on covariates?

Yes, that's what I meant when I said you could also model the
transitions as a multi-logit. There are examples in the manual on
how to code these. I don't know if people use a one-hot state vector
as a predictor or just write K different regressions.

Either way, the forward algorithm will work for general state-space
models.

- Bob

[kent]

Ben,

 

Just out of curiosity, by “the elk study”, do you mean Morales et al. Ecology 85:2436; 2004?

 

Kent

--

[Bob Carpenter]

Yes, it's the data from Morales et al. (2004). I just followed the instructions
for getting data from the vignette for the R moveHMM (Michelot, Langrock,
and Patterson's package):

trackData <- read.table("http://www.esapubs.org/archive/ecol/E085/072/elk_data.txt",
sep="\t", header=TRUE)[1:735, c(1, 2, 3, 7)]

But I also simulated data. With only four elks, this isn't going to
be much use as a hierarchical model example. Are there larger data
sets I could use, maybe with richer sets of predictors (like maybe something
categorical per animal like sex I can use in a hierarchical model)?
Ideally, with longer than 175 steps---that few data points in a mixture
doesn't exactly lead to tight posteriors.

- Bob

[kent]

There probably are people with more data, but we only had data from four elk when we wrote the paper. I haven’t been involved in analyzing data like these since then. I work on evolutionary biology and population genetics.

Kent

[Ben Lambert]

Kent - apologies for referring to it as the 'elk' paper! A very good and useful read.

There's a collaborator of mine who may have data (hourly distance and angle data) for 20 elephants. I need to discuss it with her, but she may be open to allowing the data to be made public (post publication, that is.)

Bob - thanks again for your efforts here. I have now managed to code the 'backwards' algorithm, and all seems to be working ok. The next stage is the introduction of covariates for switching probabilities.

Best,

Ben 

[Bob Carpenter]

Great---it's harder going backward because of the index
fiddling, and you have to save the forward elements!

Would you mind sharing the code?

A simple way to check the forward-backward algorithm is that
the column log_sum_exp values should all be the total likelihood
(up to arithmetic precision).

For two states, it'd all be more efficient if it were hard
coded with bernoulli transitions and a (0,1) constrained variable
for each state rather than a pair of simplexes with discrete
transitions.

- Bob

> On Jul 6, 2016, at 7:21 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Kent - apologies for referring to it as the 'elk' paper! A very good and useful read.
>

> There's a student of mine who may have data (hourly distance and angle data) for 20 elephants. I need to discuss it with her, but she may be open to allowing the data to be made public (post publication, that is.)

>
> Bob - thanks again for your efforts here. I have now managed to code the 'backwards' algorithm, and all seems to be working ok. The next stage is the introduction of covariates for switching probabilities.
>
> Best,
>
> Ben
>

[Kent Holsinger]

No need to  apologize. I was flattered to learn that it's "canonical."

Kent Holsinger

Professor of Biology

On Jul 6, 2016, at 7:21 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:

Kent - apologies for referring to it as the 'elk' paper! A very good and useful read.

There's a student of mine who may have data (hourly distance and angle data) for 20 elephants. I need to discuss it with her, but she may be open to allowing the data to be made public (post publication, that is.)

Bob - thanks again for your efforts here. I have now managed to code the 'backwards' algorithm, and all seems to be working ok. The next stage is the introduction of covariates for switching probabilities.

Best,

Ben 

--

[Ben Lambert]

Hi Bob,

The full code is below (I basically got the backwards part from the Stan manual, and - to be honest - I don't really understand it! Any clarification on how it works would be great.) Using it on fake data, it seems to do a pretty good job of returning the state vector. 

Also, I want to record the log-likelihood for each data point so I can calculate WAIC etc. Do you know how I can do this?

Best,

Ben

functions{

  ## Function that returns the log pdf of the wrapped-Cauchy

  real wrappedCauchy(real aPhi, real aRho, real aMu) {

    return(- log(2*pi()) + log((1-aRho^2)/(1+aRho^2-2*aRho*cos(aPhi-aMu))));

  }

  

  ## Function that returns a random sample from wrapped-Cauchy by first

  ## sampling from the Cauchy, then using modular arithmetic

  real wrappedCauchy_rng(real aRho, real aMu){

    real aCauchyRand;

    real aWCRand;

    

    aCauchyRand = cauchy_rng(pi()+aMu,-log(aRho));

    if (aCauchyRand>0){

      aWCRand = fmod(aCauchyRand,2*pi()) - pi();

    } else{

      aWCRand = -fmod(-aCauchyRand,2*pi()) + pi();   

    }

    return(aWCRand);

  }

  

}

data {

  int<lower=0> N;

  vector<lower=-pi(), upper=pi()>[N] phi;

  vector<lower=0>[N] l;

  int<lower=1> K;

}

parameters {

  simplex[K] theta[K];

  vector<lower=-pi(), upper=pi()>[K] mu;

  real<lower=0,upper=1> rho[K];

  positive_ordered[K] alpha;

  positive_ordered[K] beta;

}  

model {

  vector[K] log_theta[K];

  vector[K] lp;

  for (k in 1:K) log_theta[k] = log(theta[k]);

  lp = rep_vector(-log(K), K); ## Initial state log-probability

  for (n in 1:N) {

    for (k in 1:K){

      lp[k]

        = log_sum_exp(log_theta[k] + lp)

        + wrappedCauchy(phi[n], rho[k], mu[k])

        + weibull_lpdf(l[n]| alpha[k], beta[k]);

    }

  }

  target += log_sum_exp(lp);

  alpha[1] ~ normal(5,5);

  alpha[2] ~ normal(20,5);

  beta[1] ~ normal(200,5);

  beta[2] ~ normal(1000,5);

  rho[1] ~ beta(2,10);

  rho[2] ~ beta(8,10);

  

}

generated quantities {

  int<lower=1,upper=K> y_star[N];

  real log_p_y_star;

  

  ## Do backwards algorithm to get the estimated states

  {

    int back_ptr[N, K];

    real best_logp[N, K];

    real best_total_logp;

    for (k in 1:K)

      best_logp[1, K] = wrappedCauchy(phi[1], rho[k], mu[k])

                      + weibull_lpdf(l[1] | alpha[k], beta[k]);

    for (t in 2:N) {

      for (k in 1:K) {

      best_logp[t, k] = negative_infinity();

        for (j in 1:K) {

          real logp;

          logp = best_logp[t-1, j]

          + log(theta[j, k]) + wrappedCauchy(phi[t], rho[k], mu[k])

                             + weibull_lpdf(l[t] | alpha[k], beta[k]);

          if (logp > best_logp[t, k]) {

            back_ptr[t, k] = j;

            best_logp[t, k] = logp;

          }

        }

      }

    }

    log_p_y_star = max(best_logp[N]);

    for (k in 1:K)

      if (best_logp[N, k] == log_p_y_star)

      y_star[N] = k;

      for (t in 1:(N - 1))

      y_star[N - t] = back_ptr[N - t + 1,

      y_star[N - t + 1]];

  }

}

[Ben Lambert]

Also, in the attached paper the authors state that we can get better estimates of the state vector, by `rao-blackwellizing'; taking the average of the probabilities of each state, rather than the states themselves. Is this possible for the above code? 

Best,

Ben

[Bob Carpenter]

Yes, the log likelihood is just the bit you're incrementing
target += with! The rest is the prior.

You should work through how the forward and backward algorithm
work if you're going to continue in this area. It's not that
hard and is a good exercise in marginalization!

- Bob

[Bob Carpenter]

That's what's being done implicitly. If you coded it up
directly in JAGS with discrete variables for the states,
that'd be the un-Rao-Blackwellized version. (It's basically
just a matter of whether you marginalize out or sample
some variable.)

- Bob

[Ben Lambert]

Ok, thanks. 

I know that the line: 

target += log_sum_exp(lp);

Increments the log probability, but this is not pointwise (is it?) In which case, I don't know how to get the log probability per each data point, and thus how to calculate WAIC etc. Or am I way off?

Best,

Ben

[Bob Carpenter]

Right, that's the total log density. I'm not
sure how do to pointwise log density because it's
a time series, not i.i.d. observations. What exactly
do you want to calculate?

- Bob

[Ben Lambert]

Hi Bob,

I was wanting to calculate WAIC or estimate LOO-CV (using the `loo' package.) As you say, am not sure how to go about this due to its time series nature. Anyone with any ideas here?

An alternative way to proceed would be to manually partition the data into a training set (for the first T observations) and a test set (for the last t < T) observations. This may be the best thing to do here to be honest, albeit at a time cost.

Best,

Ben

[Bob Carpenter]

People often do prediction one time point at a time in these
cases, but that is very intensive.

I believe Aki, Andrew, or Jonah have already commented on this,
but if you want more feedback, I'd suggest a new thread with
a more LOO/WAIC-oriented title.

- Bob

[Ben Lambert]

Done. Thanks again Bob - you really saved our bacon here! I will mark this question 'complete'.