Correlated Random Walk with Measurement Error: Trouble Converging

Modeling
7

You will certainly want to restrict your angle parameter

vector<lower=-pi(), upper=pi()>[T - 1] turn_angle;

Otherwise the posterior will be highly multi-modal as turn_angle could just be shifted by 2*pi. Furthermore, putting a uniform prior on an unconstrained parameter is asking for trouble.

The disadvantage of restricting the range between -pi and pi is that the sampler cannot easily move across the boundary. Yet, as the parameter actually describes an angle this should be easily possible on a circular space.

Another approach would be to use a unit vector which is then transformed to polar coordinates (shortly discussed in section 19.4 of the Stan manual). I have often found it more effective to use an unconstrained 2-d vector instead with a soft constraint on its length, i.e. len ~ normal(1, 0.1);

The following code illustrates this approach for your model:

functions{
  //Function that returns the log pdf of the wrapped-Cauchy
  real wrappedCauchy_lpdf(real y, real rho, real mu) {
    return(log(1/(2*pi())) + log((1-rho^2)/(1+rho^2-2*rho*cos(y-mu))));
    }
}

data {
  int<lower=0> T;   // length of the time series
  vector[T] xo;
  vector[T] yo;
}

parameters {
  real<lower=0, upper=1> rho;
  real<lower=-pi(), upper=pi()> m;
  real<lower=0> shape;
  real<lower=0> scale;
  real<lower=0> sigma;
  vector<lower=0>[T - 1] steps;
  // unconstrained 2-D vectors
  vector[2] angle_raw[T - 1];
}

transformed parameters{
  vector[T] x;
  vector[T] y;
  vector<lower=-pi(), upper=pi()>[T-1] turn_angle;
  vector<lower=0>[T-1] angle_len;

  // transform to polar coordinates
  for (t in 2:T) {
    turn_angle[t - 1] = atan2(angle_raw[t - 1, 1], angle_raw[t - 1, 2]);
    angle_len[t - 1] = sqrt(angle_raw[t - 1, 1]^2 + angle_raw[t - 1, 2]^2);
  }
  
  x[1] = 0;
  y[1] = 0; 
  for(t in 2:T){
    x[t] = x[t - 1] + cos(turn_angle[t - 1]) * steps[t - 1];
    y[t] = y[t - 1] + sin(turn_angle[t - 1]) * steps[t - 1];
    }
}

model {
  // soft constraint on vector length
  angle_len ~ normal(1, 1e-1);
  // Jacobian correction for polar coordinates
  target += - sum(log(angle_len));

  turn_angle[1] ~ uniform(-pi(), pi());
  m ~ normal(0,0.5);
  rho ~ beta(4,1);
  scale ~ normal(0,1);
  shape ~ gamma(2,1);
  sigma ~ normal(0, .25);
  
  for(t in 2:T){
    steps[t - 1] ~ weibull(shape, scale);
    if (t > 2)
      turn_angle[t - 1] ~ wrappedCauchy(rho, turn_angle[t - 2] + m);
  }
    
  for (t in 1:T) {
    xo[t] ~ normal(x[t], sigma);
    yo[t] ~ normal(y[t], sigma);
  }

}

Note the Jacobian correction to account for the change to polar coordinates. For me this parameterization gets rid of divergent trajectories and low BMI warnings.

1 Like