Gaussian process in a plane



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

Cov_exp_quad on a 2D grid: rstan error

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.



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.