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.