Skip to content

Instantly share code, notes, and snippets.

@VictorGarritano
Created November 27, 2017 11:53
Show Gist options
  • Save VictorGarritano/129628b520166927df74158204a337a6 to your computer and use it in GitHub Desktop.
Save VictorGarritano/129628b520166927df74158204a337a6 to your computer and use it in GitHub Desktop.
import numpy as np
import matplotlib.pyplot as plt
n = 50
Xtest = np.linspace(-5, 5, n).reshape(-1,1)
def kernel(a, b, param):
sqdist = np.sum(a**2, 1).reshape(-1,1) + np.sum(b**2, 1) - 2*np.dot(a, b.T)
return np.exp(-0.5 * (1/param) * sqdist)
param = 0.1
K_ss = kernel(Xtest, Xtest, param)
L = np.linalg.cholesky(K_ss + 1e-15*np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n,3)))
# plt.plot(Xtest, f_prior)
# plt.title('Three samples from the GP prior')
# plt.grid()
# plt.xlim([-5, 5])
# plt.show()
Xtrain = np.array([-4, -3, -2, -1, 1]).reshape(5, 1)
ytrain = np.sin(Xtrain)
K = kernel(Xtrain, Xtrain, param)
L = np.linalg.cholesky(K + 0.00005*np.eye(len(Xtrain)))
K_s = kernel(Xtrain, Xtest, param)
Lk = np.linalg.solve(L, K_s)
mu = np.dot(Lk.T, np.linalg.solve(L, ytrain)).reshape((n,))
s2 = np.diag(K_ss) - np.sum(Lk**2, axis=0)
stdv = np.sqrt(s2)
L = np.linalg.cholesky(K_ss + 1e-6*np.eye(n) - np.dot(Lk.T, Lk))
f_post = mu.reshape(-1,1) + np.dot(L, np.random.normal(size=(n,3)))
plt.plot(Xtrain, ytrain, 'bs', ms=8)
plt.plot(Xtest, f_post)
plt.gca().fill_between(Xtest.flat, mu-2*stdv, mu+2*stdv, color="#dddddd")
plt.plot(Xtest, mu, 'r--', lw=2)
# pl.axis([-5, 5, -3, 3])
plt.title('Three samples from the GP posterior')
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment