apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-02 14:45:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge remote-tracking branch 'Falkor/newGP' into newGP

Browse source
This commit is contained in:
Ricardo Andrade 2013-01-29 16:49:44 +00:00
commit 81917ceba5
4 changed files with 47 additions and 85 deletions
Download patch file Download diff file Expand all files Collapse all files

10
GPy/examples/poisson.py
View file

@ -25,16 +25,14 @@ seed=default_seed
"""
X = np.arange(0,100,5)[:,None]
F = np.round(np.sin(X/18.) + .1*X)
E = np.random.randint(-3,3,20)[:,None]
F = np.round(np.sin(X/18.) + .1*X) + np.arange(5,25)[:,None]
E = np.random.randint(-5,5,20)[:,None]
Y = F + E
pb.plot(X,F,'k-')
pb.plot(X,Y,'ro')
pb.figure()
likelihood = GPy.inference.likelihoods.poisson(Y,scale=6.)
likelihood = GPy.inference.likelihoods.poisson(Y,scale=1.)
m = GPy.models.GP(X,likelihood=likelihood)
#m = GPy.models.GP(data['X'],Y=likelihood.Y)
#m = GPy.models.GP(X,Y=likelihood.Y)
m.constrain_positive('var')
m.constrain_positive('len')

27
GPy/inference/EP.py
View file

@ -48,13 +48,13 @@ class EP:
self.tau_tilde = np.zeros(self.N)
self.v_tilde = np.zeros(self.N)
def restart_EP(self):
"""
Set the EP approximation to initial state
"""
self.tau_tilde = np.zeros(self.N)
self.v_tilde = np.zeros(self.N)
self.mu = np.zeros(self.N)
def _compute_GP_variables(self):
#Variables to be called from GP
mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model
sigma_sum = 1./self.tau_ + 1./self.tau_tilde
mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2
Z_ep = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant
return self.tau_tilde[:,None], mu_tilde[:,None], Z_ep
class Full(EP):
def fit_EP(self):
@ -122,12 +122,7 @@ class Full(EP):
self.np1.append(self.tau_tilde.copy())
self.np2.append(self.v_tilde.copy())
#Variables to be called from GP
mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model
sigma_sum = 1./self.tau_ + 1./self.tau_tilde
mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2
Z_ep = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant
return self.tau_tilde[:,None], mu_tilde[:,None], Z_ep
return self._compute_GP_variables()
class DTC(EP):
def fit_EP(self):
@ -212,7 +207,8 @@ class DTC(EP):
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N
self.np1.append(self.tau_tilde.copy())
self.np2.append(self.v_tilde.copy())
return self.tau_tilde[:,None], self.v_tilde[:,None], self.Z_hat[:,None], self.tau_[:,None], self.v_[:,None]
return self._compute_GP_variables()
class FITC(EP):
def fit_EP(self):
@ -313,4 +309,5 @@ class FITC(EP):
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N
self.np1.append(self.tau_tilde.copy())
self.np2.append(self.v_tilde.copy())
return self.tau_tilde[:,None], self.v_tilde[:,None], self.Z_hat[:,None], self.tau_[:,None], self.v_[:,None]
return self._compute_GP_variables()

84
GPy/inference/likelihoods.py
View file

@ -9,65 +9,18 @@ import pylab as pb
from ..util.plot import gpplot
class likelihood:
def __init__(self,Y,location=0,scale=1):
"""
Likelihood class for doing Expectation propagation
"""
Likelihood class for doing Expectation propagation
:param Y: observed output (Nx1 numpy.darray)
..Note:: Y values allowed depend on the likelihood used
"""
:param Y: observed output (Nx1 numpy.darray)
..Note:: Y values allowed depend on the likelihood used
"""
def __init__(self,Y,location=0,scale=1):
self.Y = Y
self.N = self.Y.shape[0]
self.location = location
self.scale = scale
def plot1D(self,X,mean,var,Z=None,mean_Z=None,var_Z=None,samples=0):
"""
Plot the predictive distribution of the GP model for 1-dimensional inputs
:param X: The points at which to make a prediction
:param Mean: mean values at X
:param Var: variance values at X
:param Z: Set of points to be highlighted in the plot, i.e. inducing points
:param mean_Z: mean values at Z
:param var_Z: variance values at Z
:samples: Number of samples to plot
"""
assert X.shape[1] == 1, 'Number of dimensions must be 1'
gpplot(X,mean,var.flatten())
if samples: #NOTE why don't we put samples as a parameter of gpplot
s = np.random.multivariate_normal(mean.flatten(),np.diag(var),samples)
pb.plot(X.flatten(),s.T, alpha = 0.4, c='#3465a4', linewidth = 0.8)
#pb.subplot(211)
#self.plot1Da(X,mean,var,Z,mean_Z,var_Z)
def plot1Da(self,X,mean,var,Z=None,mean_Z=None,var_Z=None):
"""
Plot the predictive distribution of the GP model for 1-dimensional inputs
:param X_new: The points at which to make a prediction
:param Mean_new: mean values at X_new
:param Var_new: variance values at X_new
:param X_u: input (inducing) points used to train the model
:param Mean_u: mean values at X_u
:param Var_new: variance values at X_u
"""
assert X.shape[1] == 1, 'Number of dimensions must be 1'
gpplot(X,mean,var.flatten())
pb.errorbar(Z.flatten(),mean_Z.flatten(),2*np.sqrt(var_Z.flatten()),fmt='r+')
pb.plot(Z,mean_Z,'ro')
"""
def plot1Db(self,X_obs,X,phi,Z=None):
assert X_obs.shape[1] == 1, 'Number of dimensions must be 1'
gpplot(X,phi,np.zeros(X.shape[0]))
pb.plot(X_obs,(self.Y+1)/2,'kx',mew=1.5)
pb.ylim(-0.2,1.2)
if Z is not None:
pb.plot(Z,Z*0+.5,'r|',mew=1.5,markersize=12)
"""
def plot2D(self,X,X_new,F_new,U=None):
"""
Predictive distribution of the fitted GP model for 2-dimensional inputs
@ -106,6 +59,10 @@ class probit(likelihood):
L(x) = \\Phi (Y_i*f_i)
$$
"""
def __init__(self,Y,location=0,scale=1):
assert np.sum(np.abs(Y)-1) == 0, "Output values must be either -1 or 1"
likelihood.__init__(self,Y,location,scale)
def moments_match(self,i,tau_i,v_i):
"""
Moments match of the marginal approximation in EP algorithm
@ -146,6 +103,10 @@ class poisson(likelihood):
L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
$$
"""
def __init__(self,Y,location=0,scale=1):
assert len(Y[Y<0]) == 0, "Output cannot have negative values"
likelihood.__init__(self,Y,location,scale)
def moments_match(self,i,tau_i,v_i):
"""
Moments match of the marginal approximation in EP algorithm
@ -203,20 +164,19 @@ class poisson(likelihood):
sigma2_hat = m2 - mu_hat**2 # Second central moment
return float(Z_hat), float(mu_hat), float(sigma2_hat)
def plot1Db(self,X,X_new,F_new,F2_new=None,U=None):
pb.subplot(212)
#gpplot(X_new,F_new,np.sqrt(F2_new))
pb.plot(X_new,F_new)#,np.sqrt(F2_new)) #FIXME
pb.plot(X,self.Y,'kx',mew=1.5)
if U is not None:
pb.plot(U,np.ones(U.shape[0])*self.Y.min()*.8,'r|',mew=1.5,markersize=12)
def predictive_mean(self,mu,variance):
return np.exp(mu*self.scale + self.location)
def predictive_variance(self,mu,variance):
return mu
def _log_likelihood_gradients():
raise NotImplementedError
def plot(self,X,phi,X_obs,Z=None):
assert X_obs.shape[1] == 1, 'Number of dimensions must be 1'
gpplot(X,phi,np.zeros(X.shape[0]))
pb.plot(X_obs,self.Y,'kx',mew=1.5)
if Z is not None:
pb.plot(Z,Z*0+pb.ylim()[0],'k|',mew=1.5,markersize=12)
class gaussian(likelihood):
"""
Gaussian likelihood

11
GPy/models/sparse_GP.py
View file

@ -35,9 +35,13 @@ class sparse_GP(GP):
:type beta: float
:param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
:type normalize_(X|Y): bool
:parm likelihood: a GPy likelihood, defaults to gaussian
:param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1
:param powerep: power-EP parameters [$\eta$,$\delta$], defaults to [1.,1.]
:type powerep: list
"""
def __init__(self,X,Y=None,kernel=None, X_uncertainty=None, beta=100., Z=None,Zslices=None,M=10,normalize_X=False,normalize_Y=False,likelihood=None,method_ep='DTC',epsilon_ep=1e-3,epsilon_em=.1,power_ep=[1.,1.]):
def __init__(self,X,Y=None,kernel=None,X_uncertainty=None,beta=100.,Z=None,Zslices=None,M=10,normalize_X=False,normalize_Y=False,likelihood=None,method_ep='DTC',epsilon_ep=1e-3,power_ep=[1.,1.]):
if Z is None:
self.Z = np.random.permutation(X.copy())[:M]
@ -53,7 +57,7 @@ class sparse_GP(GP):
self.has_uncertain_inputs=True
self.X_uncertainty = X_uncertainty
GP.__init__(self, X=X, Y=Y, kernel=kernel, normalize_X=normalize_X, normalize_Y=normalize_Y,likelihood=likelihood,epsilon_ep=epsilon_ep,epsilon_em=epsilon_em,power_ep=power_ep)
GP.__init__(self, X=X, Y=Y, kernel=kernel, normalize_X=normalize_X, normalize_Y=normalize_Y,likelihood=likelihood,epsilon_ep=epsilon_ep,power_ep=power_ep)
self.trYYT = np.sum(np.square(self.Y)) if not self.EP else None
@ -91,6 +95,9 @@ class sparse_GP(GP):
else:
if self.hetero_noise:
print "rick's stuff here"
else:
self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices).sum()
self.psi1 = self.kern.K(self.Z,self.X)