Difference between Stan estimates and those from moveHMM for an animal movement model

[Ben Lambert]

Hi,

I have been doing some work with animal movement models, and Bob kindly helped me out with coding this up in Stan (see code at bottom of this message.) I mentioned before in passing that I was seeing some differences in the estimates of the latent state from the R package 'moveHMM' and from my Stan model when using 'optimizing' to find the posterior modes. Since I am setting uniform priors on parameters, I figured that the results should be the same.

In this message I provide a working example to highlight the different results from each approach using fake data. (Was hoping this might help someone else in the future.)

Below I use the moveHMM package to generate the fake data (with an exponential step length distribution and a wrapped Cauchy angle distribution), then estimate the model,

library(moveHMM)

stepPar <- c(1,10) # mean1, mean2, sd1, sd2

anglePar <- c(0,0,0.3,0.7) # mean1, mean2, k1, k2

stepDist <- "exp"

angleDist <- "wrpcauchy"

beta <- t(c(-3,-3))

lTest <- simData(nbAnimals=1,nbStates=2,stepDist=stepDist,angleDist=angleDist,stepPar=stepPar,

                anglePar=anglePar,obsPerAnimal=1000,beta=beta,states = TRUE)

plot(lTest)

## Fit model, plot results and extract states

aFit <- fitHMM(data=lTest,nbStates = 2,stepDist=stepDist,angleDist=angleDist,stepPar=stepPar,

               anglePar=anglePar)

plot(aFit)

lStates <- viterbi(aFit)

## Determine correlation with true states

cor(lTest$states[-c(1,1000)],lStates[-c(1,1000)])

The correlation between the true state and the fitted is pretty good, usually above 0.95. If I do the same with my Stan code, I get different latent states. Note below I remove the first and last observations since Stan doesn't like NAs,

lTest_short <- lTest[-c(1,1000),]

fit <- optimizing(aModel,data=list(l=lTest_short$step,K=2,N=dim(lTest_short)[1],phi=lTest_short$angle))

lStates_MAP <- fit$par[11:1008]

And if I look at the correlation between these two it varies from run to run, but it barely gets above a magnitude of 0.6. The Stan model tends to get stuck in a state, meaning that it is not good at reproducing reality when there are rapid changes in the latent state.

Does anyone know why there might be this difference between the two approaches? Thanks!

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))));

  }

}

data {

  int<lower=0> N;

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

  vector<lower=0>[N] l;

  int<lower=1> K;

}

parameters {

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

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

  positive_ordered[K] alpha;

  simplex[K] theta[K];

}  

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])

        + exponential_lpdf(l[n] | alpha[k]);

  }

  target += log_sum_exp(lp);

}

generated quantities {

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

  real log_p_y_star;

  {

    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])

                      + exponential_lpdf(l[1] | alpha[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])

                      + exponential_lpdf(l[t] | alpha[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 Goodrich]

On Saturday, September 3, 2016 at 7:40:21 PM UTC-4, Ben Lambert wrote:

Does anyone know why there might be this difference between the two approaches? Thanks!

I do not, but one thing to try is to see if you get the same answer using optim() with method = "L-BFGS-B". Just create a stanfit object called sf by calling stan() with iter = 0 and then call optim() with

• par = runif(get_num_upars(sf), min = -2, max = 2)
  
• fn = log_prob
  
• gr = grad_log_prob
• ... as adjust_transform = FALSE and object = sf
• method = "L-BFGS-B"
• control = list(fnscale = -1)

That will separate differences between your Stan program and moveHMM from differences in how optimization works in Stan and R. There is nothing at all special about the optimizers in Stan except that they can utilize autodiff.

Ben

[Bob Carpenter]

In the model:

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])
+ exponential_lpdf(l[n] | alpha[k]);
}
target += log_sum_exp(lp);

the part highlighted with ^^^^ isn't correct.

Thinking in terms of a N x K matrix for lp, what you want in
general for an HMM is this:

lp[n + 1, k] = SUM_k' lp[n, k'] // previous state density
+ log_theta[k', k] // transition probability
+ log p_emission(y[n] | theta(k)) // emission density

...

target += log_sum_exp(lp[N]); // final column for last time N

You can do this with two lp vectors in a rolling fashion, but you
can't do it with just one.

Your y[n] = (phi[n], l[n]) and your theta[k] = (rho[k], mu[k], alpha[k])
and the p_emission

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

[Ben Lambert]

Hi Bob,

Thanks for your message, I think I understand what you mean by comparison with the forwards part of the algorithm on page 3 here: http://www.cs.columbia.edu/~mcollins/fb.pdf 

So I have changed my code to the following,

model {

  matrix[K,K] log_theta;

  matrix[N,K] lp;

  for (kk in 1:K){

    for (k in 1:K){

      log_theta[kk,k] = log(theta[kk][k]);

    }

  }

    

  lp[1] = to_row_vector(rep_vector(-log(K), K)); ## Initial state log-probability

  

  for (n in 2:N) {

    for(kk in 1:K){

      real aTempLP;

      aTempLP = 0;

      for (k in 1:K){

        aTempLP

          = aTempLP + lp[n-1,k] + log_theta[kk,k];

      }

      lp[n,kk] = aTempLP + wrappedCauchy(phi[n], rho[kk], mu[kk])

          + exponential_lpdf(l[n] | alpha[kk]);

    }

  }

  target += log_sum_exp(lp[N]);

}

However I now get the following warning which freezes the optimisation,

Exception thrown at line 57: stan::math::exponential_log: Inverse scale parameter is 0, but must be > 0!

This is the like where I call exponential_lpdf. Any idea where I'm going wrong? I have spent the past hour debugging by print, but can't seem to find the source.

Many thanks,

Ben

[Bob Carpenter]

What's happening is that alpha[kk] is becoming zero at some
point. Are you sure the model has a mode with the
data you provide? Sometimes the likelihood can grow toward
a bound and cause underflow.

I'd suggest putting a prior on alpha[kk] that bounds it away
from zero like a lognormal().

Or you could try sampling rather than optimizing.

There are examples in the manual, too. You only need the
forward part of forward backward for Stan. If you use
all of forward-backward, you can often do analytical updates
of theta and weighted updates of the emission densities,
which would be faster, but is harder to code.

Also, your vectorization for log_theta was fine as it was
before but it'll be just as efficient with loops.

- Bob

[Ben Lambert]

Placing lognormal priors on alpha doesn't seem to alleviate the problem.

I'm pretty sure something is awry with the model as it currently stands. moveHMM seems able to find a maximum ok, so assume that the mode is findable. (Also the latent states estimated from the model are almost 100% in agreement with the true states.)

Even with sampling there are a huge number of complaints about alpha being near zero. When I print the log probability lp[N] I get huge values: -10^300!

It looks like lp[i] is basically doubling (in magnitude) each step because of the sum over lp[i-1]. Does this sound sensible?

> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

This was my fault for giving you a buggy version of the
model in the first place. Sorry about that. I got too
excited after the conference in Seattle. Then didn't
pay enough attention to answering your email query.

Here's a model that works with von-Mises turns and Weibull
steps. You'll see that the theta matrix needs to be transposed,
then log_sum_exp-ed, because the forward algorithm adds all
the ways you can get from previous state to next state, so
that requires log_sum_exp (and the transposition, so the
simplex constraint is applied in the right place).


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_tr[K];
vector[K] lp;
vector[K] lp_p1;

for (k_from in 1:K)

for (k in 1:K)

log_theta_tr[k, k_from] = log(theta[k_from, k]);
lp = rep_vector(-log(K), K);

for (n in 1:N) {

for (k in 1:K)

lp_p1[k]
= log_sum_exp(log_theta_tr[k] + lp)
+ von_mises_lpdf(turn[n] | mu[k], kappa[k])
+ weibull_lpdf(dist[n] | alpha[k], sigma[k]);
lp = lp_p1;
}
target += log_sum_exp(lp);

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

For the two-state HMM, I simulated

theta <- matrix(c(0.90, 0.10,
0.10, 0.90), 2, 2, byrow=TRUE);

mu <- c(-0.2, 0)
kappa <- c(0.5, 4)

alpha <- c(1, 3);
sigma <- c(0.25, 0.5);

and gave it a single 10,000-step simulated path. It converged
and fit fine in 200 iterations total in about 4 minutes and
recovered:

Inference for Stan model: move-hmm.
1 chains, each with iter=200; warmup=100; thin=1;
post-warmup draws per chain=100, total post-warmup draws=100.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1,1] 0.90 0.00 0.00 0.89 0.90 0.90 0.90 0.91 40 1.03
theta[1,2] 0.10 0.00 0.00 0.09 0.10 0.10 0.10 0.11 40 1.03
theta[2,1] 0.11 0.00 0.00 0.10 0.10 0.11 0.11 0.12 55 1.01
theta[2,2] 0.89 0.00 0.00 0.88 0.89 0.89 0.90 0.90 55 1.01
mu[1] -0.25 0.00 0.05 -0.34 -0.28 -0.24 -0.21 -0.16 100 0.99
mu[2] -0.01 0.00 0.01 -0.02 -0.02 -0.01 -0.01 0.00 100 0.99
kappa[1] 0.45 0.00 0.02 0.41 0.43 0.44 0.46 0.49 42 0.99
kappa[2] 3.95 0.01 0.09 3.79 3.89 3.97 4.01 4.11 100 1.00
alpha[1] 1.00 0.00 0.01 0.97 0.99 1.00 1.01 1.02 100 1.00
alpha[2] 3.04 0.00 0.04 2.97 3.01 3.04 3.06 3.11 87 1.01
sigma[1] 0.25 0.00 0.00 0.24 0.25 0.25 0.25 0.26 80 0.99
sigma[2] 0.50 0.00 0.00 0.50 0.50 0.50 0.50 0.51 91 0.99
lp__ -11312.40 0.28 2.09 -11317.05 -11313.46 -11312.18 -11310.94 -11309.50 57 1.01

Samples were drawn using NUTS(diag_e) at Sun Sep 4 21:53:50 2016.


So I think this one's OK. As you can see, the state
transitions are well estimated, but the locations of
the von Mises have fairly broad posteriors, which makes
sense given their wide scales.

It would be even better to fit this model with proper
circular parameters. I'll need to figure out how to
do that with our unit_vector parameter type, but I
should be able to get that done for a case study (a
unit_vector[2] will provide a point on a circle which
can then be translated back to an angle---all this trig
makes my head hurt).

When I get some time write the case study out, I'll also
show how to pull out the best predicted path of states
(use the Viterbi algorithm) and the marginals over states
at each point in time (that requires full forward-backward).
Great you found Michael Collins's notes---he's a great
teacher.

Thanks for your patience.

- Bob

[Ben Lambert]

Hi Bob,

Please don't apologise -- I am very grateful to the time you personally have given me on this forum (I know I take more than most.) Also thanks very much for the new code -- looks pretty good!

I have now implemented the forwards-backwards algorithm (I think) and seem to be getting reasonable performance. The correlation between the estimated state and the true state is about 0.94. This is a bit lower than the moveHMM package estimates (a correlation of 0.97), but not far off.

Suspect there may be a small error in the code somewhere, but this isn't doing too badly. I have changed the step dist to an exponential since this is more in line with my data, but otherwise your code remains intact.

If anyone gets the chance to look through the below and comment on whether I'm making a silly mistake somewhere, then let me know. 

Best,

Ben

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;

}

transformed parameters {

  vector[K] log_theta_tr[K];

  

   for (k_from in 1:K)

    for (k in 1:K)

      log_theta_tr[k, k_from] = log(theta[k_from, k]);

}

model {

  vector[K] lp;

  vector[K] lp_p1;

  

  lp = rep_vector(-log(K), K);

  

  # Forwards algorithm

  for (n in 1:N) {

    for (k in 1:K){

      lp_p1[k]

        = log_sum_exp(log_theta_tr[k] + lp)

        + von_mises_lpdf(turn[n] | mu[k], kappa[k])

        + exponential_lpdf(dist[n] | alpha[k]);

    }

    lp = lp_p1;

  }

  target += log_sum_exp(lp);

}

generated quantities{

  matrix[N,K] B;

  vector[K] lp_pTotal;

  vector[K] lp_p2;

  vector[K] lp_p1a;

  vector[K] lp_a;

  matrix[N,K] A;

  matrix[N,K] qState;

  vector[N] pState;

  

  # Forwards algorithm

  lp_a = rep_vector(-log(K), K);

  

  for (n in 1:N) {

    for (k in 1:K){

      lp_p1a[k]

        = log_sum_exp(log_theta_tr[k] + lp_a)

        + von_mises_lpdf(turn[n] | mu[k], kappa[k])

        + exponential_lpdf(dist[n] | alpha[k]);

    }

    lp_a = lp_p1a;

    A[n] = to_row_vector(lp_a);

  }

  

  ## Backwards algorithm http://www.cs.columbia.edu/~mcollins/fb.pdf

  B[N] = rep_row_vector(1,K);

  lp_pTotal = rep_vector(0,K);

  

  for(n in 1:(N-1)){

   for(k in 1:K){

     lp_p2[k]

        = log_sum_exp(log_theta_tr[k] + lp_pTotal)

        + von_mises_lpdf(turn[N - n + 1] | mu[k], kappa[k])

        + exponential_lpdf(dist[N - n + 1] | alpha[k]);

   } 

   lp_pTotal = lp_p2;

   B[N - n] = to_row_vector(lp_pTotal);

  }

  

  ## Estimate unnormalised state probs

  for(n in 1:N){

    for(k in 1:K){

      qState[n,k] = A[n,k] + B[n,k];

    }

  }

 

  # Normalise these to find prob of state 2 (am using a 2 state model)

  for(n in 1:N){

      pState[n] = exp(qState[n,2] - log_sum_exp(qState[n]));

  }

}

[Bob Carpenter]

Assuming the models are the same, you should get exactly
the same fit. I assume moveHMM also marginalizes out the
discrete states to do optimization. But I also recall them
doing something to compute the stationary distribution for
initialization rather than taking a uniform distribution.

I also take it you're comparing the best-path states computed
via the Viterbi algorithm.

I can take a look at the rest of the code when I get some
more time.

What I'd do is simulate data and see how it fits in both
systems.

- Bob

[Ben Lambert]

I have made a look at the moveHMM package -- they use the viterbi algorithm to compute the most likely sequence of states, rather than forwards-backwards (see https://cran.r-project.org/web/packages/moveHMM/vignettes/moveHMM-guide.pdf). 

This is the right thing to do, it seems. The forwards-backwards can compute the most likely individual state for any given point in time, whereas the viterbi is what we want to use if we want to compute the most likely sequence of states (see https://en.wikipedia.org/wiki/Forward%E2%80%93backward_algorithm). The differences between these two is quite subtle!

When I replace the forwards-backwards code (in the 'generated quantities' section) with the viterbi I get the same performance as moveHMM in estimating the state probabilities! Looks like all is working now, thanks again,

generated quantities {

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

  real log_p_y_star;

  {

    int back_ptr[N, K];

    real best_logp[N, K];

    real best_total_logp;

    for (k in 1:K)

      best_logp[1, K] = von_mises_lpdf(turn[1] | mu[k], kappa[k])

                      + exponential_lpdf(dist[1] | alpha[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]) + von_mises_lpdf(turn[t] | mu[k], kappa[k])

                      + exponential_lpdf(dist[t] | alpha[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]];

  }

}

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

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

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Ben Lambert]

I also found this to be a handy reference on the differences between viterbi and the forwards-backwards,

https://www.robots.ox.ac.uk/~vgg/rg/slides/hmm.pdf

[Bob Carpenter]

Glad to hear the performance is the same.

The Viterbi algorithm finds the most probable state sequence z*
given the observed data y (emissions) and parameters theta
(transition and distributions of emissions per state):

z* = ARGMAX_z p(z | y, theta)

whereas forward-backward computes the marginal probability
that state z[t] = k, given the observed data and parameters

Pr[z[t] = k | y, theta]

You can use that to find the most likely state at time t,
but that's not necessarily the same as z*[t].

With forward backward and a little more work, you can compute
the probabitity of any subsequence of states. And with Viterbi,
you can generalize to output the N-best possible sequences in
order.

- Bob

> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --

> 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.

[Bob Carpenter]

Larry's tutorial's a classic---I saw him present it more
than once at Bell Labs when I used to work in natural language
and speech recognition. Bishop's machine learning book is
nice and it contains the generalized algorithms and approximations,
as is Durbin's book on biological sequences. I imagine
the Langrock et al. HMM book is good---I like their doc for
moveHMM.

- Bob

> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
>
> --

> 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.

[Ben Lambert]

Yes, am now doing cross-validation to estimate the out-of-sample predictive capability of a few variants of the models. Following Aki's advice in another thread, am going to try to estimate the probability of the data for blocks,

1      2     3     4  |  5

1      2     3  |  4  |  5

1      2  |  3  |  4     5

1   |  2  |  3     4     5

1   |  2     3     4     5

where ``|'' indicates a border between a training and testing set. (I guess this is a sort of K-Fold cross validation.)

Am using the forwards and backwards algorithms to determine the predictive probability of the test sequences. Will let you know how I get on!

Best,

Ben

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

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

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

You want to use just the forward algorithm starting
from the end of the observed sequence to evaluate
the predictive likelihood of the test data.

I didn't get what the ones with two bars mean.

- Bob

> > > 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.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > To unsubscribe from this group and all its topics, 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.
> >
> >
> >
> > --
> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.

[Ben Lambert]

Yes, sorry I mis-typed. Meant just the forward.

> > > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> > > To post to this group, send email to stan-...@googlegroups.com.
> > > For more options, visit https://groups.google.com/d/optout.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

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

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

And it's not the full forward, but the forward *after* all
the observations you see. So if you're predicting state 4--5,
it's just the conditional likelihood of those states you want
to compute.

- Bob

> > > > 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.
> > >
> > > --
> > > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > > To unsubscribe from this group and all its topics, 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.
> > >
> > >
> > >
> > > --
> > > 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.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > To unsubscribe from this group and all its topics, 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.
> >
> >
> > --
> > 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.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.

[Ben Lambert]

Yes, I have been using the state probabilities for the last training set observation to initialise my forward algorithm on the test set. (Rather than in the full forward where we assign equal probability to each state initially.) 

This means that I use the forwards-backwards algorithm on the training set to determine the state probabilities for the last observation. I then use these to initialise a forward run on the test set. This means that I do the following and use ``lp_testOverall'' to measure the out-of-sample predictive capability of the model,

generated quantities{

  matrix[N,K] B;

  vector[K] lp_pTotal;

  vector[K] lp_p2;

  vector[K] lp_p1a;

  vector[K] lp_a;

  matrix[N,K] A;

  matrix[N,K] qState;

  vector[N] pState;

  matrix[N,K] pStateAll;

  

  ## Test sample

  vector[K] lp_test;

  vector[K] lp_testTemp;

  real lp_testOverall;

  

  # Forwards algorithm

  lp_a = rep_vector(-log(K), K);

  

  for (n in 1:N) {

    for (k in 1:K){

      lp_p1a[k]

        = log_sum_exp(log_theta_tr[k] + lp_a)

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

        + exponential_lpdf(dist[n] | alpha[k]);

    }

    lp_a = lp_p1a;

    A[n] = to_row_vector(lp_a);

  }

  

  

  ## Backwards algorithm http://www.cs.columbia.edu/~mcollins/fb.pdf

  B[N] = rep_row_vector(1,K);

  lp_pTotal = rep_vector(0,K);

  

  for(n in 1:(N-1)){

   for(k in 1:K){

     lp_p2[k]

        = log_sum_exp(log_theta_tr[k] + lp_pTotal)

        + wrappedCauchy(turn[N - n + 1] , rho[k], mu[k])

        + exponential_lpdf(dist[N - n + 1] | alpha[k]);

   } 

   lp_pTotal = lp_p2;

   B[N - n] = to_row_vector(lp_pTotal);

  }

  

  ## Estimate unnormalised state probs

  for(n in 1:N){

    for(k in 1:K){

      qState[n,k] = A[n,k] + B[n,k];

    }

  }

  

  # Normalise these to find prob of state 2 (am using a 2 state model)

  for(n in 1:N){

      pState[n] = exp(qState[n,2] - log_sum_exp(qState[n]));

      pStateAll[n] = exp(qState[n] - log_sum_exp(qState[n]));

  }

  

  # Out-of-sample predictive capability estimation

  ## Initialise state probs according to last training set obs

  lp_test = to_vector(pStateAll[N]);

 ## Use forward algorithm 

  for (n in 1:N1) {

    for (k in 1:K){

      lp_testTemp[k]

        = log_sum_exp(log_theta_tr[k] + lp_test)

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

        + exponential_lpdf(distTest[n] | alpha[k]);

    }

    lp_test = lp_testTemp;

  }

  ### sum over X_k to marginalise

  lp_testOverall = log_sum_exp(lp_test);

  

}

> > > > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> > > > To post to this group, send email to stan-...@googlegroups.com.
> > > > For more options, visit https://groups.google.com/d/optout.
> > >
> > > --
> > > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > > To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

> > > To post to this group, send email to stan-...@googlegroups.com.
> > > For more options, visit https://groups.google.com/d/optout.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> > To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

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

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

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/72BEgxLIZjo/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

> On Sep 7, 2016, at 10:56 AM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Yes, I have been using the state probabilities for the last training set observation to initialise my forward algorithm on the test set. (Rather than in the full forward where we assign equal probability to each state initially.)
>
> This means that I use the forwards-backwards algorithm on the training set to determine the state probabilities for the last observation. I then use these to initialise a forward run on the test set. This means that I do the following and use ``lp_testOverall'' to measure the out-of-sample predictive capability of the model,

The key here is that you only want to include the predicted
data likelihoods. I think you can do this with only the
forward algorithm. But I don't have time to write it all out
now.

- Bob

[Ben Lambert]

Yes that is what I am doing.

[Ben Lambert]

Just in case someone will find this useful for future work, I have now extended the model to allow for covariates to affect the probability of switching states. Specifically, I have used the same framework explained in https://cran.r-project.org/web/packages/moveHMM/vignettes/moveHMM-guide.pdf, on pages 19 and 20, where the off-diagonal probabilities are multinomial logits.

The below seems to reproduce the parameters that I use to generate the data, and be in line with the states estimated via fitHMM itself. Note: I have checked the model only for the 2- and 3-state models, with a single covariate. 

Does anyone think that it'd be worthwhile added functionality to Stan that implements estimation of these animal movement HMMs? That is, rather than the user having to code them up as I've done below. As a starter it might be good to try to add the distributions and models from the moveHMM package. Could always extend to more general HMMs later.

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] turn;

  vector<lower=0>[N] dist;

  int<lower=1> K;

  

  ## Out-of-sample test

  int N1;

  vector<lower=-pi(), upper=pi()>[N1] turnTest;

  vector<lower=0>[N1] distTest;

  

  # Covariate

  int nCovs;

  matrix[N,nCovs+1] X;

}

parameters {

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

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

  positive_ordered[K] alpha;

  matrix[K*(K-1),nCovs+1] beta;

}

transformed parameters {

  ## Gamma here is a matrix of the log probabilities not probabilities as in the moveHMM doc

  matrix[K,K] Gamma[N];

  matrix[K,K] Gamma_tr[N];

  

  for(n in 1:N){

    int aCount;

    aCount = 1;

    for(k_from in 1:K){

      for(k in 1:K){

        if(k_from==k){

          Gamma[n,k_from,k] = 1;

      }else{

          Gamma[n,k_from,k] = exp(beta[aCount] * to_vector(X[n]));

          aCount = aCount + 1;

      }

    }

    }

    Gamma[n] = log(Gamma[n]/sum(Gamma[n]));

  }

  

  for(n in 1:N)

    for(k_from in 1:K)

      for(k  in 1:K)

        Gamma_tr[n, k, k_from] = Gamma[n,k_from,k];

  

}

model {

  vector[K] lp;

  vector[K] lp_p1;

  

  lp = rep_vector(-log(K), K);

  

  # Forwards algorithm

  for (n in 1:N) {

    for (k in 1:K){

      lp_p1[k]

        = log_sum_exp(to_vector(Gamma_tr[n,k]) + lp)

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

        + exponential_lpdf(dist[n] | alpha[k]);

    }

    lp = lp_p1;

  }

  target += log_sum_exp(lp);

}

generated quantities{

## Viterbi to estimate most likely state vector

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

  real log_p_y_star;

  {

    int back_ptr[N, K];

    real best_logp[N, K];

    real best_total_logp;

    for (k in 1:K)

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

        + exponential_lpdf(dist[1] | alpha[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]

            + Gamma[t,j,k] + wrappedCauchy(turn[t] , rho[k], mu[k])

            + exponential_lpdf(dist[t] | alpha[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)

      state[N] = k;

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

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

      state[N - t + 1]];

  }

}

[Ben Goodrich]

On Saturday, September 10, 2016 at 12:07:52 PM UTC-4, Ben Lambert wrote:

Does anyone think that it'd be worthwhile added functionality to Stan that implements estimation of these animal movement HMMs? That is, rather than the user having to code them up as I've done below. As a starter it might be good to try to add the distributions and models from the moveHMM package. Could always extend to more general HMMs later.

 
I think whenever you are targeting the functionality of an existing R package that is well-known in a field, then it makes sense to build a package using rstan::rstan_package_skeleton(), which will do it like rstanarm does. Or ask the moveHMM maintainer to do it so that Stan is an estimation option.

Ben

[Bob Carpenter]

There are some separate issues here.

1. Add new distributions, like

real wrappedCauchy_lpdf(real phi | reals rho, reals mu)

2. Add new general primitive functions

vector change_bearing(vec x, vec y);

3. Add some HMM-specific densities

real hmm(vector[] y | vector[] theta, vector[] phi, function Foo);

The arguments must be arrays of vectors making up
the observation sequence y, transition simplexes theta,
and state-specific parameters phi, where Foo is going to
have to be a function that generates the observation densities
for state k using parameters phi[k].

4. Write a high-level R package like moveHMM that wraps up whole models
some options. For that, you want to look at RStanArm and the R package
skeleton to which Ben was referring. If you can prespecify a finite
number of Stan models (perhaps conditional) following RStanArm, it
makes the derived package much simpler to install than RStan itself.

- Bob


> On Sep 10, 2016, at 12:07

> --

> 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.