Constraints fail on transformed parameter block in hierarchical HMM model

Modeling
,
1

Hi All,
I am trying to expand Hidden Markov Models to fit to typical movement data. The starting point is based on earlier discussions between @bgoodri and @Ben_Lambert. One extension that I am trying is to fit one hierarchical model to several animals, allowing the parameters that govern movement (ie step length and turning angle) to vary between the individuals.
The issue that I am running into is that initialization fail despite setting constraints to parameters in the transformed parameter block.
Here is the model I tried: hmm_covariates_hier.stan (5.1 KB)

And some code to generate data:

library(moveHMM)
library(rstan)

n_obs <- 200
ss <- simData(nbAnimals = 5,nbStates = 2,stepDist = "exp",
              angleDist = "wrpcauchy",stepPar = c(2,0.1),anglePar = c(0,1,0.15,0.75),
              obsPerAnimal = n_obs)
data <- list(N=(n_obs-2)*5,turn=ss$angle[!is.na(ss$angle)],
         dist=ss$step[!is.na(ss$angle)],K=2,nCovs=1,
         X=matrix(rep(1,(n_obs-2)*5),ncol=1),S=5,Ind_ID=ss$ID[!is.na(ss$angle)])

m <- stan("hmm_covariates_hier.stan",data=data,control=list(adapt_delta=0.99))
#this result in many initial error:


Error evaluating the log probability at the initial value.
validate transformed params: rho_ind[3] is -4.38754, but must be greater than or equal to 0
Rejecting initial value:
Error evaluating the log probability at the initial value.
validate transformed params: mu_ind[5] is -4.66398, but must be greater than or equal to -3.14159

Initialization between (-2, 2) failed after 100 attempts. 
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

I tried to play a bit around setting an initialization function, but this did not solve anything.
Am I missing something?

2

So far all that’s missing is the error message.

3

Thanks for your fast answer, updated the original post.

4

So I can’t open a .Stan file on my phone easily but that’s telling you that you have a (transformed?) parameter with a constraint that is not being met. Are you missing an explanation transformation in an assignment?

5

Ugh, exp not explanation

6

It looks like both rho_ind, alpha_ind, and mu_ind are all multiplied by unbounded parameters (error), which would cause them to take values outside of the constraints specified by their declarations:

parameters {
  ...
  matrix[S,K] error[3];
  ...
}

transformed parameters {
  ...
  matrix<lower=-pi(),upper=pi()>[S,K] mu_ind;
  matrix<lower=0>[S,K] rho_ind;
  matrix<lower=0>[S,K] alpha_ind;
  ...
  mu_ind = rep_matrix(mu + rep_vector(sigma_mu, K), S)' .* error[1,,];
  rho_ind = rep_matrix(rho + rep_vector(sigma_rho, K), S)' .* error[2,,];
  alpha_ind = rep_matrix(alpha + rep_vector(sigma_alpha, K), S)' .* error[3,,];
  ..
}
model {
  ...
  to_vector(error[1,,]) ~ normal(0,1);
  to_vector(error[2,,]) ~ normal(0,1);
  to_vector(error[3,,]) ~ normal(0,1);
  ...
}

So, as @sakrejda suggests, maybe using some kind of link function or transformation to satisfy the lower and upper bounds of support would be a good solution.

7

Ah yes, good point. Here is how I re-wrote the model:

transformed parameters {
...
  mu_ind = sigma_mu * error[1,,]; //mu_ind is N(0,sigma_mu)
  mu_ind = (2 * pi() / (max(mu_ind) - min(mu_ind))) * (mu_ind - max(mu_ind)) + pi(); //re-scaling between   -pi and pi
  mu_ind = mu_ind + rep_matrix(mu,S)'; //mu_ind is N(mu,sigma_mu) and scaled between -pi ad pi

  rho_ind = rep_matrix(rho + rep_vector(sigma_rho, K), S)' .* exp(error[2,,]);

  alpha_ind = rep_matrix(alpha + rep_vector(sigma_alpha, K), S)' .* exp(error[3,,]);
...
}

The model runs despite some initialization failures for the mu_ind parameter, I guess that I am somehow doing the math wrong there.

Now the issue is that the chains show all kind of bad behaviors (no mixing, non-stationarity for some of them, divergence, issue with treedepth …). So I wonder if this is an issue with this parametrization, or if this particular type of model is tricky and I should feed in more (simulated) data.

8

If it’s not a supervised HMM, you run into all the same problems as with other mixture models. You might want to take a look at Michael Betancourt’s mixture case study.

Initialization failure winds up happening either due to missing constraints (some values should have non-zero density but don’t) or due to numerical issues.