# Multivariate Gaussian Process Regression Inference

Hi,
I’m new to Gaussian Processes and I want to extend the model from chapter 17.4 to the multivariate case.

Here is what I came up with:

``````model_string_gp <-"
data {
int<lower=1> D;
int<lower=1> N1;
int<lower=1> N2;
vector[D] x1[N1];
vector[N1] y1;
vector[D] x2[N2];
}
transformed data {
int<lower=1> N;
vector[D] x[N1+N2];
vector[N1+N2] mu;
matrix[N1+N2, N1+N2] Sigma;
//cov_matrix[N1+N2] Sigma;
N = N1 + N2;
for (n in 1:N1) x[n] = x1[n];
for (n in 1:N2) x[N1 + n] = x2[n];
for (i in 1:N) mu[i] = 0;
for (i in 1:N)
for (j in 1:N)
Sigma[i, j] = exp(-dot_self(x[i] - x[j]))
+ i == j ? 0.1 : 0.0;
}
parameters {
vector[N2] y2;
}
model {
vector[N] y;
for (n in 1:N1) y[n] = y1[n];
for (n in 1:N2) y[N1 + n] = y2[n];
y ~ multi_normal(mu, Sigma);
}
"
``````

I changed “cov_matrix[N1+N2] Sigma” to a simple matrix because I was getting errors being not symmetrical.
When I fit this model to some data, I run into these errors:

``````Rejecting initial value:
Error evaluating the log probability at the initial value.

Rejecting initial value:
Error evaluating the log probability at the initial value.

Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler\$call_sampler(args_list[[i]]) : Initialization failed."
``````

I changed the “init_r” but it changed nothing. The univariate model from the docs runs fine.

Hey Mick

Thanks for noting the positive definite error otherwise I woulda been lost on this one haha. I changed:

``````Sigma[i, j] = exp(-dot_self(x[i] - x[j]))
+ i == j ? 0.1 : 0.0;
``````

to

`Sigma[i, j] = exp(-dot_self(x[i] - x[j])) + (i == j ? 0.1 : 0.0);`

And the model went to working for me (notice the parenthesis around the diagonal term). I think this might be an error in the manual.

(paging @Bob_Carpenter) For the code

`````` Sigma[i, j] = exp(-dot_self(x[i] - x[j]))
+ i == j ? 0.1 : 0.0;
``````

The `+` gets evaluated before the `==` right (this code is straight from page 250 of the 2.15 manual)? I guess this is all gonna be replaced by @rtrangucci 's stuff in the next manual?

Mick,

are you talking about Gaussian random fields?

Andre

thanks, that seems to fix the problem. I changed the parenthesis and it worked. I guess that needs to be fixed in the docs as I just copy/pasted this line from there.

I don’t know what you mean, I just extended the example from chapter 17.4 on Gaussian Processes to the multivariate case.

Yes, @rtrangucci rewrote the entire chapter on GPs for the next version of the manual, so at least any old problems will go away.

This is not what you want:

``````Sigma[i, j] = exp(-dot_self(x[i] - x[j])) + i == j ? 0.1 : 0.0;
``````

because the `?` and `:` bind too losely. To code this, you want to set

``````Sigma[i, j] = exp(-dot_self(x[i] - x[j])) + i == j ? 0.1 : 0.0;
``````

in the loop and then use a separate loop to increment the diagonal as

``````for (k in 1:K)
Sigma[k, k] = Sigma[k, k] + 0.1;
``````

In the next release, we’ll have compound add/assign: `+=`

thanks Bob, do you mean

`Sigma[i, j] = exp(-dot_self(x[i] - x[j]))`

in the loop and then

``````for (k in 1:K)
Sigma[k, k] = Sigma[k, k] + 0.1;
``````

in another loop?

What is the difference to the parenthesis suggestion from @bbbales2 ?

The second apparently leads to less bugs :D

Also, for radial basis function Gaussian processes, check out the function `cov_exp_quad` in the docs. It lets you build a covariance matrix with one line of code. It’s really handy cause sometimes code for GPs can get out of hand (there’s some examples using it here: https://github.com/stan-dev/stancon_talks/tree/master/2017/Contributed-Talks/08_trangucci).

In that case, if you want to add stuff to the diagonal, you’ll need to use the second loop as well. `cov_exp_quad` might be a bit faster too.

Fewer comparisons and less complicated expressions.