Gaussian process in a plane

Modeling
1

This is a really simple question. I have a set of x and y coordinates as data. I perform some manipulations on them which are defined in a transformed parameters block, and then in the model block I feed them to cov_exp_quad. After trying several different data structures the following code snippet works:

transformed parameters {
  vector[N] xc = x - x_c;
  vector[N] yc = y - y_c;
  vector[N] yhat = -xc * sin(phi) + yc * cos(phi);
  vector[N] xhat = -(xc * cos(phi) + yc * sin(phi))/ci;
  row_vector[2] xyhat[N];
  xyhat[:, 1] = to_array_1d(xhat);
  xyhat[:, 2] = to_array_1d(yhat);
}

Except for some model specific details the rest just follows the examples in Michael Betancourt’s and Rob Trangucci’s tutorials.

The question is whether this data structure is as efficient as possible, and if not what is? It works in the sense that cov_exp_quad is correctly calculating the Euclidean distance between points (x,y) and (x’, y’).

2

I’m afraid there’s not a more convenient or efficient way to go between all these different structures.

You can save some time by not recaclulating sin(phi) and cos(phi). If phi is a scalar, you’re also better off with

vector[N] xhat;
{
  real cos_phi = cos(phi);  // these are all local variables
  real sin_phi = sin(phi);
  real neg_inv_ci = -inv(ci);
  xhat = cos_phi * neg_inv_ci * xc + sin_phi * neg_inv_ci * yc;
}

But the basic point is to not duplicate computation.

There may be some clever trig thing that simplifies the calculation of xhat and yhat together I haven’t thought about. That’d be the real win—if calculating yhat makes it easier to calcuate xhat or vice-versa.

3

Thanks.

phi is a scalar parameter. So is ci. The manipulations are doing a rotation and a deprojection, and there are certainly other ways to code the transformations. Also, a 2D unit vector instead of cos(phi) and sin(phi) should work. I realized recently that my first effort in that direction failed because I didn’t give the parameter an informative prior, which is really necessary. I need to experiment some more.

By the way this is an example of a model with parameters defined on a circle. If I’m ever satisfied with it maybe I’ll write up a case study.

4

We’d love to get a real example of circular stats. I’ve been meaning to code up a proper movement HMM case study with circular Cauchy and a unit vector type, but haven’t gotten around to it yet.

5

I’ve been playing around with various ways to code my model and finally came up with this: there’s a rotate and a stretch (I think that’s the term of art) in there, and somewhat to my surprise it’s substantially faster (10-40% in a limited number of trials so far) to write the transformations as matrix multiplies. Here’s my current transformed parameters block:

transformed parameters {
  real sin_i = sqrt(1.-cos_i^2);
  vector[N] xc = x - x_c;
  vector[N] yc = y - y_c;
  real sin_phi = sin(phi);
  real cos_phi = cos(phi);
  matrix[N, 2] X;
  matrix[N, 2] Xhat;
  matrix[2, 2] Rot;
  matrix[2, 2] Stretch;
  vector[N] r;
  vector[N] xhat;
  vector[N] yhat;
  row_vector[2] xyhat[N];
  
  X = append_col(xc, yc);
  Rot = [ [-cos_phi, -sin_phi],
          [-sin_phi,  cos_phi]];
  Stretch = diag_matrix([1./cos_i, 1.]');
  Xhat = X * Rot * Stretch;
  xhat = Xhat[ : , 1];
  yhat = Xhat[ : , 2];
  r = sqrt(xhat .* xhat + yhat .* yhat);
  
  xyhat[:, 1] = to_array_1d(xhat);
  xyhat[:, 2] = to_array_1d(yhat);
}

phi and cos_i are scalar parameters. Even though phi is an angle I found it better not to make it bounded. I use a principled informative prior instead. Since it enters the model through sin(phi), cos(phi) pairs I could use a 2D unit vector instead. That turns out to be slower. cos_i is constrained to (0,1). Why not make the parameter sec_i which would be bounded below by 1? That’s slower too. cos_i also needs an informative prior.

We’d love to get a real example of circular stats.

I’ll try to write something up, but it may be a while. In the meantime I’ll get code and some sample data on github before I have to drop this for a while.

Thanks for the suggestions.

1 Like