Skip to content

Instantly share code, notes, and snippets.

@osdf
Created May 19, 2011 06:30
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save osdf/980292 to your computer and use it in GitHub Desktop.
Interactive sequential linear regression
#!/usr/bin/env python
"""
__author__: Martin Felder
Demo for sequential bayesin learning
for linear regression models. See
Chapter 3.3 of 'Pattern recognition
and machine learning' by Chris Bishop, 2001
"""
import scipy as sp;
import numpy as npy;
from numpy import mat
from pylab import *
from copy import copy
def bivar_gauss(u, v, m, Sinv):
x = array([u,v])
return exp(-0.5*dot(dot((x-m).transpose(),Sinv),x-m))
def mesh2vector(U, V):
"""convert given mesh to a vector with v[:,i] giving gridpoint i's coords"""
return mat([U.flatten(),V.flatten()])
def vector2mesh(vec, template):
"""convert given mesh to a vector with v[:,i] giving gridpoint i's coords"""
return vec.reshape(size(template,axis=0),size(template,axis=1))
def multivar_gauss(X, m, Sinv):
"""computes a multivariate Gaussian for multiple sample vectors (cols of X)"""
nVec = size(X,axis=1)
expo = empty(nVec)
for i in range(nVec):
z = X[:,i]-m
expo[i] = z.getT()*Sinv*z
# normalize
g = exp(-0.5*expo)
return g/g.sum()
class basisFunctionAffine:
def getNoisyVal(self,*args):
"""0=a0, 1=a1, 2,3=random range, 4=#points, 5=noise"""
nPoints = args[4]
x = npy.random.uniform(args[2],args[3],nPoints)
y = self.dotProd(x,args[0],args[1]) + npy.random.normal(0,0.2,nPoints)
return x, y
def dotProd(self,x,*args):
"""returns dot product of basis function with given vector"""
return args[0] + args[1]*x
def designMatrix(self,X):
"""returns design matrix of basis function"""
return npy.mat([ones(len(X)),X]).transpose()
class bayesLinReg:
def __init__(self, alpha, beta):
"""inititalize with a Gaussian prior and no measurements"""
self.alpha = alpha
self.beta = beta
self.initPrior()
self.setBasisFct('basisFunctionAffine')
def initPrior(self):
self.pMean = mat(zeros([2,1]))
self.pCov = mat(identity(2) / self.alpha)
def prior(self,U,V):
if self.pCov[0,1] != self.pCov[1,0]:
print "Error: Asymmetric prior covariance!"
sys.exit(1)
return vector2mesh(multivar_gauss(mesh2vector(U,V), self.pMean, self.pCov.getI()), U)
#bivariate_normal(U, V, sigmax=sqrt(self.pCov[0,0]), sigmay=sqrt(self.pCov[1,1]), mux=self.pMean[0], muy=self.pMean[1], sigmaxy=sqrt(self.pCov[0,1]))
def likelihood(self,x,t,U,V):
"""likelihood function for a single observation"""
#K = len(U)
#lk = npy.empty_like(V)
std = 1.0/self.beta
lk = normpdf(t, self.basisFct.dotProd(x,U,V), std)
return lk
def setBasisFct(self, fctname):
self.basisFct = eval(fctname+'()')
def getPosteriorParams(self,X,T):
"""return posterior mean and inverse(!) covariance depending on samples given (Bishop eqn. 3.50/51)"""
phi = self.basisFct.designMatrix(X)
SnI = self.pCov.getI() + self.beta*(phi.getT()*phi)
mn = SnI.getI()*(self.pCov.getI()*self.pMean + self.beta*phi.getT()*mat(T).getT())
return mn, SnI
def posterior(self,X,T,U,V):
"""return posterior depending on samples given (Bishop eqn. 3.49)"""
mn, SnI = self.getPosteriorParams(X,T)
return vector2mesh(multivar_gauss(mesh2vector(U,V), mn, SnI), U)
#return npy.nan_to_num(bivariate_normal(U, V, sigmax=sqrt(Sn[0,0]), sigmay=sqrt(Sn[1,1]), mux=mn[0], muy=mn[1], sigmaxy=sqrt(Sn[0,1])))
#return normpdf(U,mn[0],Sn[0,0])*normpdf(V,mn[1],Sn[1,1])
def drawW(self,X,T):
"""create a design matrix from given samples"""
mn, SnI = self.getPosteriorParams(X,T)
return npy.random.multivariate_normal(array(mn).flatten(), SnI.getI())
noise = 0.2 # stddev of noise in the data
beta = (1.0/noise)**2 # precision parameter of data
alpha = 2.0 # precision of prior
trueW = array([-0.3,0.5])
# generate the data
f = basisFunctionAffine()
X = []
Y = copy(X)
x1 = array([-1,1]) # for plotting
# prepare the prior/posterior density
du, dv = 0.02, 0.02
u = arange(-1.0, 1.0, du)
v = arange(-1.0, 1.0, dv)
U,V = meshgrid(u, v)
# setup Bayesian regression
bayes = bayesLinReg(alpha, beta)
bayes.setBasisFct('basisFunctionAffine')
# setup graphics
fig = figure(figsize=[18,5])
pLike = fig.add_subplot(131)
pLike.set_title('likelihood', fontsize=20)
pPost = fig.add_subplot(132)
pPost.set_title('prior density', fontsize=20)
pPost.set_xlabel('$w_0$', fontsize=15)
pPost.set_ylabel('$w_1$', fontsize=15)
pDat = fig.add_subplot(133)
pDat.set_title(r'data space', fontsize=20)
nP = 0
# plot the prior density
p = bayes.prior(U, V)
pPost.imshow(p, interpolation='bilinear', origin='lower',\
extent=(-1,1,-1,1))
# plot true data
pDat.plot(x1, f.dotProd(x1,trueW[0],trueW[1]), color='g',alpha=1.0,linewidth=3)
pDat.hold(True)
# plot the data
for i in range(6):
w1 = npy.random.uniform(-1,1)
w2 = npy.random.uniform(-1,1)
pDat.plot(x1, f.dotProd(x1,w1,w2), color='r',alpha=0.5,linewidth=3)
pDat.set_ylim(-1,1)
pDat.set_xlabel('$x$', fontsize=15)
pDat.set_ylabel('$t$', fontsize=15)
pDat.hold(False)
class InteractiveBayes(object):
def onclick(self, event):
if not event.button==1: return True
# add the click locations
X.append(event.xdata)
Y.append(event.ydata)
nP = len(X)
print nP, event.xdata, event.ydata
# plot the likelihood
p = bayes.likelihood(X[nP-1],Y[nP-1],U,V)
im = pLike.imshow(p, interpolation='bilinear', origin='lower',\
extent=(-1,1,-1,1))
pLike.hold(True)
pLike.plot([trueW[0]],[trueW[1]],'+',markersize=10,markeredgecolor='w',markeredgewidth=2)
pLike.set_xlabel('$w_0$', fontsize=15)
pLike.set_ylabel('$w_1$', fontsize=15)
hold(False)
# plot the posterior density
p = bayes.posterior(X[0:nP],Y[0:nP], U, V)
im = pPost.imshow(p, interpolation='bilinear', origin='lower',\
extent=(-1,1,-1,1))
pPost.hold(True)
pPost.plot([trueW[0]],[trueW[1]],'+',markersize=10,markeredgecolor='w',markeredgewidth=2)
pPost.set_xlabel('$w_0$', fontsize=15)
pPost.set_ylabel('$w_1$', fontsize=15)
pPost.set_title('posterior density', fontsize=20)
hold(False)
# plot true data
pDat.plot(x1, f.dotProd(x1,trueW[0],trueW[1]), color='g',alpha=1.0,linewidth=3)
pDat.hold(True)
# plot the data
for i in range(6):
w1, w2 = bayes.drawW(X[0:nP],Y[0:nP])
pDat.plot(x1, f.dotProd(x1,w1,w2), color='r',alpha=0.5,linewidth=3)
pDat.plot(X[0:nP],Y[0:nP],'o',markeredgecolor='b', markerfacecolor='w',markeredgewidth=1.5)
#plot(u,p.sum(axis=0)*10-1,linewidth=2,color='g') # plot marginalized gaussian
pDat.set_ylim(-1,1)
pDat.set_title(r'data space', fontsize=20)
pDat.set_xlabel('$x$')
pDat.set_ylabel('$t$')
pDat.hold(False)
fig.canvas.draw()
bay = InteractiveBayes()
fig.canvas.mpl_connect('button_press_event', bay.onclick)
show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment