# Simulate response from 2d GP

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;
real<lower=0> alpha;
real<lower=0> sigma;
}

transformed data {
matrix[N, N] cov1 =   cov_exp_quad(x[,1], alpha, rho)
+ 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, rho)
+ 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.

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

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

You could rescale the input vectors.

``````transformed data {
vector x_r[N];
for (i in 1:N) {
x_r[i,1] = x[i,1] / rho;
x_r[i,2] = x[i,2] / rho;
}
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, rho)
+ diag_matrix(rep_vector(1e-10, N));

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

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