example_draw.py

from sklearn.gaussian_process.kernels import RBF

# First we set up the training data
xs_train = training_points # The points that will reign in our prediction
ys_train = train_temp # The corresponding temperatures for those points

# Then the space for which we are trying to predict 
xs_test = np.atleast_2d(np.linspace(0, 20, 1000)).T

# We will use the SKLearn RBF function
def rbf_kernel(x, y, alpha): 
    rbf = RBF(alpha)
    return rbf.__call__(x, y)

# Calculate the covariance matrix using the test data
K_ss = rbf_kernel(xs_test, xs_test, 2.)

# Get the "standard deviation" equivalent of the covariance matrix
# Note the small constant added to make the matrix positive semi-definite
L = np.linalg.cholesky(K_ss + 1e-10 * np.eye(len(xs_test)))

# Generate some normals (5 of them) 
normals = np.random.normal(size=(len(xs_test), 5))

# Dot them with the cholesky decomp. of the covariance matrix (scale them up)
f = np.dot(L, normals)

# Plot the 5 sampled functions.
plt.plot(xs_test, f)
plt.title('Five samples from the GP prior')
plt.show()