Spatial anisotropy

Modeling
1

Good evening,

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);
  }
2

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]