Skip to content

Instantly share code, notes, and snippets.

@stober
Created Feb 16, 2013
Embed
What would you like to do?
Gaussian Process in Python
#!/usr/bin/python
"""
Author: Jeremy M. Stober
Program: GP.PY
Date: Thursday, July 17 2008
Description: Example of Gaussian Process Regression.
"""
from numpy import *
import pylab
class Kernel(object):
"""
Kernel from Bishop's Pattern Recognition and Machine Learning pg. 307 Eqn. 6.63.
"""
def __init__(self,*args):
self.thetas = args
def __call__(self,x,y):
exponential = self.thetas[0] * exp( -0.5 * self.thetas[1] * sum( (x - y)**2 ) )
linear = self.thetas[3] * dot(x,y)
constant = self.thetas[2]
return exponential + constant + linear
class OrnsteinKernel(object):
"""
Ornstein-Uhlenbeck process kernel.
"""
def __init__(self,theta):
self.theta = theta
def __call__(self,x,y):
return exp(-self.theta * sum(abs(x-y)))
def covariance(kernel, data):
return reshape([kernel(x,y) for x in data for y in data], (len(data),len(data)))
def draw_multivariate_gaussian(mean,C):
ndim = len(mean)
z = random.standard_normal(ndim)
# Better numerical stabability than cholskey decomposition for
# near-singular matrices C.
[U,S,V] = linalg.svd(C)
A = U * sqrt(S)
return mean + dot(A,z)
def train(data,kernel):
mean = zeros(len(data))
C = covariance(kernel,data)
return (mean,C)
def predict(x, data, kernel, C, t):
"""
The prediction equations are from Bishop pg 308. eqns. 6.66 and 6.67.
"""
k = [kernel(x,y) for y in data]
Cinv = linalg.inv(C)
m = dot(dot(k,Cinv),t)
sigma = kernel(x,x) - dot(dot(k,Cinv),k)
return (x,m,sigma)
#kernel = OrnsteinKernel(1.0)
kernel = Kernel(1.0, 64.0, 0.0, 0.0)
# Some sample training points.
xpts = random.rand(10) * 2 - 1
# In the context of Gaussian Processes training means simply
# constructing the kernel (or Gram) matrix.
(m,C) = train(xpts, kernel)
# Now we draw from the distribution to sample from the gaussian prior.
t = draw_multivariate_gaussian(m,C)
pylab.figure(0)
pylab.plot(xpts, t, "+")
# Instead of regressing against some known function, lets just see
# what happens when we predict based on the sampled prior. This seems
# to be what a lot of other demo code does.
# Explore the results of GP regression in the target domain.
predictions = [predict(i,xpts,kernel,C,t) for i in arange(-1,1,0.01)]
pylab.figure(1)
x = [prediction[0] for prediction in predictions]
y = [prediction[1] for prediction in predictions]
sigma = [prediction[2] for prediction in predictions]
pylab.errorbar(x,y,yerr=sigma)
pylab.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment