# Sampling from Gaussian Process Posterior

Anyone know of a Python package that both fits a Gaussian Process to data, and also lets you sample paths from the posterior? I'm interested in sampling the colorful lines on right (b) of the following picture from Rasmussen's GPML book.

The scikit-learn package gives pointwise predictions and pointwise variances, which beautifully draw the prediction and 95% confidence interval. It also gives the Cholesky decomposition, but that has dimensions equivalent to the data set (small), and I want to draw samples on a grid (larger). From sklearn:

The GPy and pyGPs packages let you optimize the hyperparameters, and make the same graph as above, but when I try to build the kernel with the optimal hyperparameter values, and then draw a multivariate normal with the covariance matrix that GPy generates (using the [kernel_object].K(X) syntax), the lines are unstable and far out of the 95% confidence interval.

Any advice? I'm just trying to get those measly wiggly lines to show up! :)

• Do you need to sample from adjusted model? Or you want just to obtain plot similar to that at figures? In the latter case calculate covariance for your covariance function and transform $\mathcal{N}(0, I_n)$ data to desired multivariate normal distribution. Feb 24, 2016 at 16:58
• Thanks, Alexey. I'd like to sample from the adjusted model. Feb 24, 2016 at 16:59
• One issue: I'm unsure exactly which kernel and which parameters are used in the posterior, after it is optimized. If I knew it, I could manually compute the covariance matrix. Feb 24, 2016 at 17:01
• Use model.posterior_samples_f function in GPy. Note, that after optimization noise variance can be significantly greater than 0, so your models by default would not be interpolation. Feb 24, 2016 at 17:34
• Or you can return posterior covariance matrix in GPy with option full_cov=true in predict function, and then able to sample from corresponding normal distribution Feb 24, 2016 at 17:49

You want to sample posterior using the data and model given.

In this case you can:

• sample from posterior normal distribution with given mean and covariance matrix - use model.predict with full_covariance=True in case;
• use built-in function model.posterior_samples_f that does the job for you.

A sample code is below:

import GPy
import numpy as np

sample_size = 5
X = np.random.uniform(0, 1., (sample_size, 1))
Y = np.sin(X) + np.random.randn(sample_size, 1)*0.05

kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)
model = GPy.models.GPRegression(X,Y,kernel, noise_var=1e-10)

testX = np.linspace(0, 1, 100).reshape(-1, 1)
posteriorTestY = model.posterior_samples_f(testX, full_cov=True, size=3)
simY, simMse = model.predict(testX)

plt.plot(testX, posteriorTestY)
plt.plot(X, Y, 'ok', markersize=10)
plt.plot(testX, simY - 3 * simMse ** 0.5, '--g')
plt.plot(testX, simY + 3 * simMse ** 0.5, '--g')

• Brilliant! Thank you, @Alexey ! Combined with model.optimize() it works just as I had hoped. Feb 24, 2016 at 20:51