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! :)