Stan Covariance Matrix

Question

I am experimenting with stan and Gaussian Processes. After some errors I found out the root of everything is a strange behavior of the function cov_exp_quad.
In particular I do not understand why it returns a matrix that is not symmetrical.
Here the stan script:

data {
  int<lower=1> N;
  real x[N];
  real y[N];
}

transformed data {
  matrix[N,N] K = cov_exp_quad(x, y, 1, 1);
  }

generated quantities {
  matrix[N, N] Cov = K;
}

and here the R code:

library(rstan)

x <- seq(0, 1, length.out = 3)

X <- expand.grid(x, x)

data_stan <- list(N = dim(X)[1],
                  x = X[, 1],
                  y = X[, 2])

fit <- stan(file = "./stan_script.stan",
            data = data_stan,
            algorithm = "Fixed_param",
            warmup = 0,
            chains = 1,
            iter = 1)

samples <- rstan::extract(fit)

drop(samples$Cov)

This is, instead the output

> drop(samples$Cov)
           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]
 [1,] 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969 0.6065307 0.6065307 0.6065307
 [2,] 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969
 [3,] 0.6065307 0.6065307 0.6065307 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000
 [4,] 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969 0.6065307 0.6065307 0.6065307
 [5,] 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969
 [6,] 0.6065307 0.6065307 0.6065307 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000
 [7,] 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969 0.6065307 0.6065307 0.6065307
 [8,] 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000 0.8824969 0.8824969 0.8824969
 [9,] 0.6065307 0.6065307 0.6065307 0.8824969 0.8824969 0.8824969 1.0000000 1.0000000 1.0000000

Why isn't this symmetrical? Thank you

score 0 · Answer 1 · answered Aug 30 '22 at 16:20

It's making a covariance matrix using the distance between the elements of each vector, not using the Euclidean distance between points in [x,y]. For example:

#Squared distance between first element of x and all elements of y
(distSq <- (data_stan$x[1] - data_stan$y)^2) 

#Function to create covariance matrix 
#Taken from https://mc-stan.org/docs/2_22/functions-reference/covariance.html
cov_exp_quad <- function(alpha,rho,dist2) (alpha^2)*exp((-1/2*(rho^2))*dist2)

#Compare to first row of covariance matrix from Stan - identical
cbind(samples$Cov[1,1,],cov_exp_quad(1,1,distSq))

Stan Covariance Matrix

1 Answers1