Additive Gaussian Processes

Hi all,

I was playing with additive Gaussian processes similar to documentation here:

but I am having issues while sampling from the conditional distribution for f1_star and f2_star.

My notebook can be found here:

with pm.Model() as model:
    sigma1 = pm.HalfCauchy("η1", beta=5)
    sigma2 = pm.HalfCauchy("η2", beta=5)
    ls1 = pm.HalfCauchy('ls1', beta=10)
    ls2 = pm.HalfCauchy('ls2', beta=10)
    cov1 = sigma1 *, ls=ls1, active_dims=[0])
    cov2 = sigma2 *, ls=ls2, active_dims=[1]) 
    gp1 =
    gp2 =
    gp = gp1 + gp2
    sigma = pm.HalfCauchy("sigma", beta=3)
    f = gp.marginal_likelihood("f", X_2D_train, y, sigma)
    trace = pm.sample(1000, chains=4)

with model:
    f_star = gp.conditional("f_star", X_2D_star)

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 1, got 2 with shape (1000, 1000).

X_2D_star has shape (1000,2). So, I think I am using additive GP wrong. Can someone please help me with this?


Hi all,

I wonder if anyone has some tips on how to do additive GP with pymc3 and sample from the conditional distributions.


What happens if you use gp.predict? Does the following work?

with model:
    pred = gp.predict(X_2D_star)

Thanks for your reply @ckrapu.

It works as long as the shape of X_2D_star is the same as the shape in X_2D_train. However, that may not be the case in the most general case. Also, the output is a list with 2 arrays of shape (1000, 1000). How do I interpret these 2 arrays?


I’m a little stuck here, though I think this thread could be very helpful. @bwengals, any idea? Here’s a minimal reproducible example

With regard to the output of gp.predict, the first output should be the predictive mean vector and the output is the covariance matrix.

1 Like

Thanks @ckrapu.

Thats how I interpreted it but was confused because both elements of the list had shape (1000, 1000) in my case. Shouldn’t the mean be a vector of (1000, 1)?


Apologies, I missed that detail. There is clearly something fishy going on with the broadcasting as the predictive mean should only have one nontrivial dimension.

Yes it should, something is fishy here. I can find some time to look into this.


Sorry about the delay. So I ran through the notebook and got it to run by specifying the cov_func argument.

gp1 =
gp2 =

So this is a bit of a gotcha, but if you check the docstring:

Init signature:
    mean_func=< object at 0x7fa7d8182130>,
    cov_func=< object at 0x7fa7d8182160>,

so the mean function was set with the cov function.