Simulate response from 2d GP

Modeling
1

Hi,

I would like to simulate a GP response from a 2d data input using different hyperparameters. As, Betancourt shows in his case study one can simulate a GP by doing

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

  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma;
}

transformed data {
  matrix[N, N] cov =   cov_exp_quad(x, alpha, rho)
                     + diag_matrix(rep_vector(1e-10, N));
  matrix[N, N] L_cov = cholesky_decompose(cov);
}

parameters {}
model {}

generated quantities {
  vector[N] f = multi_normal_cholesky_rng(rep_vector(0, N), L_cov);
  real y[N] = normal_rng(f, sigma);
}

However the above is for one input variable.

For the 2d case I was thinking something like that but I’m not sure how to combine the two L_cov (and if this makes sense).

data {
  int<lower=1> N;
  matrix x[N,2];

  real<lower=0> rho[2];
  real<lower=0> alpha[2];
  real<lower=0> sigma;
}

transformed data {
  matrix[N, N] cov1 =   cov_exp_quad(x[,1], alpha[1], rho[1])
                     + diag_matrix(rep_vector(1e-10, N));
  matrix[N, N] L_cov1 = cholesky_decompose(cov1);
  
  matrix[N, N] cov2 =   cov_exp_quad(x[,2], alpha[2], rho[2])
                     + diag_matrix(rep_vector(1e-10, N));
  matrix[N, N] L_cov2 = cholesky_decompose(cov2);
}

parameters {}
model {}

generated quantities {
  vector[N] f = multi_normal_cholesky_rng(rep_vector(0, N), L_cov);
  real y[N] = normal_rng(f, sigma);
}

In sort, how can I simulate a response y from two predictors x1 and x2?

Thanks for your time.

2

The documentation says cov_exp_quad is defined for vector input too.

(that page says the function is deprecated but that’s probably not relevant for RStan.)

I think you can use an array of vectors:

data {
  int<lower=1> N;
  vector[2] x[N];

  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma;
}

transformed data {
  matrix[N, N] cov =   cov_exp_quad(x, alpha, rho)
                     + diag_matrix(rep_vector(1e-10, N));
  matrix[N, N] L_cov = cholesky_decompose(cov);
}
3

Indeed, someone can use vector for the input variables. But how can I set different length-scale for each input variable?

4

You could rescale the input vectors.

transformed data {
  vector[2] x_r[N];
  for (i in 1:N) {
    x_r[i,1] = x[i,1] / rho[1];
    x_r[i,2] = x[i,2] / rho[2];
  }
  matrix[N, N] cov =   cov_exp_quad(x_r, alpha, 1.0)
                     + diag_matrix(rep_vector(1e-10, N));
}

Alternatively, you asked how to combine two covariance matrices. Looking at the formula on the documentation page I think element-wise multiplication should work (not matrix multiplication)

transformed data {
  matrix[N, N] cov1 =   cov_exp_quad(x[,1], alpha[1], rho[1])
                     + diag_matrix(rep_vector(1e-10, N));
  
  matrix[N, N] cov2 =   cov_exp_quad(x[,2], alpha[2], rho[2])
                     + diag_matrix(rep_vector(1e-10, N));

  matrix[N,N] cov = cov1 .* cov2;
  matrix[N,N] L_cov = cholesky_decompose(cov);
}
1 Like