Hilbert space GPs with distance matrices

Modeling
6

Hey there, sorry I don’t think I’m really following you. I’ve now written what I believe is a multivariate form of the HSGP with two dimensions, namely lon and lat in the matrix x, using the following functions:

functions {
  // eigenvalues (2D, vectorised)
  vector get_lambda(vector L, array[] int j) {
    vector[2] lambda = square((to_vector(j) * pi()) ./ (2 * L));
    return lambda;
  }
	
  // eigenvectors (2D)
  vector get_phi(vector L, array[] int j, matrix x) {
    int n_obs = rows(x);
    int D = 2;
    matrix[n_obs, D] phi;
    for (d in 1:D) {
      phi[, d] = 1 / sqrt(L[d]) * sin(j[d] * pi() * (x[, d] * L[d]) / (2 * L[d]));
    } // d
    return phi[, 1] .* phi[, 2];
  }
	
  // spectral density (2D, vectorised)
  vector get_spd(real alpha, real l, matrix omega) {
    int M = rows(omega);
    vector[M] s = square(alpha) * 2 * pi() * square(l) * exp(-0.5 * square(l) * rows_dot_self(omega));
    return s;
  }
}

Then in my model the code is this (ignore the weird indexing due to this happening inside a partial sum function:

model {
  for (i in 1:n_spp) {
  // spectral densities
  spd[i] = get_spd(exp(log_alpha[ii]), exp(log_l[ii]), sqrt(lambda));
  // for each hypervolume
  for (j in 1:n_hyp[ii]) {
  // HSGP realisations
    f[i, j][:, 1:n_obs[ii, j]] = diag_pre_multiply(exp(log_w[ii]), 
                                                   diag_post_multiply(beta[ii, j], sqrt(spd[i]))) 
                                 * phi[ii, j][:, 1:n_obs[ii, j]];
    }
  }
}

Here, the beta, spd, and phi (the latter being computed in transformed data not shown here) components follow the Riutart Mayol paper—other than beta now being a matrix instead of a vector of standard normals to account the multivariate part—and I think/hope the log_w part is the “inverse scales” as per the Stan Functions Reference for multivariate Gaussian processes. Unfortunately there’s no references there in the manual so I’m not exactly sure if I’ve done it correctly here, but I assume w is another vector of scale parameters for each of the multivariate dimensions… would love to hear @avehtari’s thoughts on whether this looks correct!

Cheers,

Matt