I am totally new to Stan (I started digging into it two days ago). I was told that Stan might help me “fit” spatial stationary anisotropic model in a more simple way and in a much quicker way. But I will admit I am struggling.

I think I got the correct way to write Sigma in the case of the exponential stationary (without anisotropy). If not please feel free to comment (and correct) any mistake. But I don’t see how to proceed from there to add the anisotropy portion/ get all the pieces.

I am hoping that someone on this forum has encountered the case where (s)he would have to a simple spatial GP with geometric anisotropy in the variance and would be willing to share his/her model write up to get me started.

Thank you very much

Sab

// the Sigma matrix
matrix[N1,N1] L;
{
matrix[N1,N1] Sigma;
for(i in 1:(N1-1)){
for(j in (i+1):N1){
Sigma[i,j] = alpha^2 * exp((-1)*phi*dist[i,j]);
Sigma[j,i] = Sigma[i,j];
}
}
for(i in 1:N1) Sigma[i,i] = sigma_sq;
for (n in 1:N1)
Sigma[n, n] = Sigma[n, n] + square(sigma);
L = cholesky_decompose(Sigma);
}

Correct me if I’m wrong, but I think anistropy just means that the spatial dimensions each have their own characteristic length scale hyperparameter. If this is the case, take a look at section 10.3 in the Stan User’s Guide. In particular, notice this model from the Automatic Relevance Determination part:

functions {
matrix L_cov_exp_quad_ARD(vector[] x,
real alpha,
vector rho,
real delta) {
int N = size(x);
matrix[N, N] K;
real sq_alpha = square(alpha);
for (i in 1:(N-1)) {
K[i, i] = sq_alpha + delta;
for (j in (i + 1):N) {
K[i, j] = sq_alpha
* exp(-0.5 * dot_self((x[i] - x[j]) ./ rho));
K[j, i] = K[i, j];
}
}
K[N, N] = sq_alpha + delta;
return cholesky_decompose(K);
}
}
data {
int<lower=1> N;
int<lower=1> D;
vector[D] x[N];
vector[N] y;
}
transformed data {
real delta = 1e-9;
}
parameters {
vector<lower=0>[D] rho;
real<lower=0> alpha;
real<lower=0> sigma;
vector[N] eta;
}
model {
vector[N] f;
{
matrix[N, N] L_K = L_cov_exp_quad_ARD(x, alpha, rho, delta);
f = L_K * eta;
}
rho ~ inv_gamma(5, 5);
alpha ~ std_normal();
sigma ~ std_normal();
eta ~ std_normal();
y ~ normal(f, sigma);
}

For a 2-D anisotropic spatial model, you would have D=2 and you may want to modify the priors for length scales, rho[1] and rho[2]