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

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.

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.

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.

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