Multivariate Gaussian Process Regression Inference

Modeling
1

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.

2

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?

3

Mick,

are you talking about Gaussian random fields?

Andre

4

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.

5

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

6

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: +=

7

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 ?

8

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.

9

Fewer comparisons and less complicated expressions.