%matplotlib inline import numpy as np import scipy as sp import matplotlib.pyplot as plt import pods # set prior variance on w alpha = 4. # set the order of the polynomial basis set degree = 5 # set the noise variance sigma2 = 0.01 def polynomial(x, degree, loc, scale): degrees = np.arange(degree+1) return ((x-loc)/scale)**degrees data = pods.datasets.olympic_marathon_men() x = data['X'] y = data['Y'] plt.plot(x, y, 'rx') scale = np.max(x) - np.min(x) loc = np.min(x) + 0.5*scale num_data = x.shape[0] num_pred_data = 100 # how many points to use for plotting predictions x_pred = np.linspace(1880, 2030, num_pred_data)[:, None] # input locations for predictions Phi_pred = polynomial(x_pred, degree=degree, loc=loc, scale=scale) Phi = polynomial(x, degree=degree, loc=loc, scale=scale) num_samples = 10 K = degree+1 for i in xrange(num_samples): z_vec = np.random.normal(size=(K, 1)) w_sample = z_vec*np.sqrt(alpha) f_sample = np.dot(Phi_pred,w_sample) plt.plot(x_pred, f_sample) K = alpha*np.dot(Phi_pred, Phi_pred.T) for i in np.arange(10): f_sample = np.random.multivariate_normal(mean=np.zeros(x_pred.size), cov=K) plt.plot(x_pred.flatten(), f_sample.flatten()) fig, ax = plt.subplots(figsize=(8,8)) im = ax.imshow(K, interpolation='none') fig.colorbar(im) K = alpha*np.dot(Phi_pred, Phi_pred.T) + sigma2*np.eye(x_pred.size) for i in xrange(10): y_sample = np.random.multivariate_normal(mean=np.zeros(x_pred.size), cov=K) plt.plot(x_pred.flatten(), y_sample.flatten()) sigma2 = 1. K = alpha*np.dot(Phi_pred, Phi_pred.T) + sigma2*np.eye(x_pred.size) for i in xrange(10): y_sample = np.random.multivariate_normal(mean=np.zeros(x_pred.size), cov=K) plt.plot(x_pred.flatten(), y_sample.flatten()) def exponentiated_quadratic(x, x_prime, variance, lengthscale): squared_distance = ((x-x_prime)**2).sum() return variance*np.exp((-0.5*squared_distance)/lengthscale**2) def compute_kernel(X, X2, kernel, **kwargs): K = np.zeros((X.shape[0], X2.shape[0])) for i in np.arange(X.shape[0]): for j in np.arange(X2.shape[0]): K[i, j] = kernel(X[i, :], X2[j, :], **kwargs) return K K = compute_kernel(x_pred, x_pred, exponentiated_quadratic, variance=1., lengthscale=10.) fig, ax = plt.subplots(figsize=(8,8)) im = ax.imshow(K, interpolation='none') fig.colorbar(im) for i in xrange(10): y_sample = np.random.multivariate_normal(mean=np.zeros(x_pred.size), cov=K) plt.plot(x_pred.flatten(), y_sample.flatten()) class GP(): def __init__(self, X, y, sigma2, kernel, **kwargs): self.K = compute_kernel(X, X, kernel, **kwargs) self.X = X self.y = y self.sigma2 = sigma2 self.kernel = kernel self.kernel_args = kwargs self.update_inverse() def update_inverse(self): # Preompute the inverse covariance and some quantities of interest ## NOTE: This is not the correct *numerical* way to compute this! It is for ease of use. self.Kinv = np.linalg.inv(self.K+self.sigma2*np.eye(self.K.shape[0])) # the log determinant of the covariance matrix. self.logdetK = np.linalg.det(self.K+self.sigma2*np.eye(self.K.shape[0])) # The matrix inner product of the inverse covariance self.Kinvy = np.dot(self.Kinv, self.y) self.yKinvy = (self.y*self.Kinvy).sum() def log_likelihood(self): # use the pre-computes to return the likelihood return -0.5*(self.K.shape[0]*np.log(2*np.pi) + self.logdetK + self.yKinvy) def objective(self): # use the pre-computes to return the objective function return -self.log_likelihood() # set covariance function parameters variance = 16.0 lengthscale = 32 # set noise variance sigma2 = 0.05 K = compute_kernel(x, x, exponentiated_quadratic, variance=variance, lengthscale=lengthscale) K_star = compute_kernel(x, x_pred, exponentiated_quadratic, variance=variance, lengthscale=lengthscale) K_starstar = compute_kernel(x_pred, x_pred, exponentiated_quadratic, variance=variance, lengthscale=lengthscale) fig, ax = plt.subplots(figsize=(8,8)) im = ax.imshow(np.vstack([np.hstack([K, K_star]), np.hstack([K_star.T, K_starstar])]), interpolation='none') # Add lines for separating training and test data ax.axvline(x.shape[0]-1, color='w') ax.axhline(x.shape[0]-1, color='w') fig.colorbar(im) def posterior_f(self, X_test): K_starstar = compute_kernel(X_test, X_test, self.kernel, **self.kernel_args) K_star = compute_kernel(self.X, X_test, self.kernel, **self.kernel_args) A = np.dot(self.Kinv, K_star) mu_f = np.dot(A.T, y) C_f = K_starstar - np.dot(A.T, K_star) return mu_f, C_f # attach the new method to class GP(): GP.posterior_f = posterior_f model = GP(x, y, sigma2, exponentiated_quadratic, variance=variance, lengthscale=lengthscale) mu_f, C_f = model.posterior_f(x_pred) fig, ax = plt.subplots(figsize=(8,8)) im = ax.imshow(C_f, interpolation='none') fig.colorbar(im) plt.plot(x, y, 'rx') plt.plot(x_pred, mu_f, 'b-') var_f = np.diag(C_f)[:, None] std_f = np.sqrt(var_f) plt.plot(x, y, 'rx') plt.plot(x_pred, mu_f, 'b-') plt.plot(x_pred, mu_f+2*std_f, 'b--') plt.plot(x_pred, mu_f-2*std_f, 'b--') from IPython.display import YouTubeVideo YouTubeVideo('ewJ3AxKclOg') def update_inverse(self): # Perform Cholesky decomposition on matrix self.R = sp.linalg.cholesky(self.K + self.sigma2*self.K.shape[0]) # compute the log determinant from Cholesky decomposition self.logdetK = 2*np.log(np.diag(self.R)).sum() # compute y^\top K^{-1}y from Cholesky factor self.Rinvy = sp.linalg.solve_triangular(self.R, self.y) self.yKinvy = (self.Rinvy**2).sum() # compute the inverse of the upper triangular Cholesky factor self.Rinv = sp.linalg.solve_triangular(self.R, np.eye(self.K.shape[0])) self.Kinv = np.dot(self.Rinv, self.Rinv.T) GP.update_inverse = update_inverse for i in xrange(10): y_sample = np.random.multivariate_normal(mean=mu_f.flatten(), cov=C_f) plt.plot(x_pred.flatten(), y_sample.flatten()) plt.plot(x, y, 'rx')