diff --git a/GPy/examples/fitc_class_01.py b/GPy/examples/fitc_class_01.py new file mode 100644 index 00000000..413d0c14 --- /dev/null +++ b/GPy/examples/fitc_class_01.py @@ -0,0 +1,53 @@ +import pylab as pb +import numpy as np +import GPy +pb.ion() +pb.close('all') + +""" +Simple 1D classification example +:param seed : seed value for data generation (default is 4). +:type seed: int +""" +seed=10000 + +data = GPy.util.datasets.toy_linear_1d_classification(seed=seed) +X = data['X'] +Y = data['Y'][:, 0:1] +Y[Y == -1] = 0 + + + +#X = np.vstack((np.random.uniform(0,10,(10,1)),np.random.uniform(7,17,(10,1)),np.random.uniform(15,25,(10,1)))) +#Y = np.vstack((np.zeros((10,1)),np.ones((10,1)),np.zeros((10,1)))) + +# Kernel object +kernel = GPy.kern.rbf(1) + GPy.kern.white(1) + +# Likelihood object +distribution = GPy.likelihoods.likelihood_functions.probit() +likelihood = GPy.likelihoods.EP(Y,distribution) + +Z = np.random.uniform(X.min(),X.max(),(10,1)) +#Z = np.array([0,20])[:,None] +print Z + +# Model definition +m = GPy.models.generalized_FITC(X,likelihood=likelihood,kernel=kernel,Z=Z,normalize_X=False) +m.set('len',2.) + +m.ensure_default_constraints() +# Optimize +#m.constrain_fixed('iip') +m.update_likelihood_approximation() +print m.checkgrad(verbose=1) +# Parameters optimization: +m.optimize() +#m.EPEM() #FIXME + +# Plot +pb.subplot(211) +m.plot_f() +pb.subplot(212) +m.plot() +print(m) diff --git a/GPy/examples/fitc_class_02.py b/GPy/examples/fitc_class_02.py new file mode 100644 index 00000000..cd8d6700 --- /dev/null +++ b/GPy/examples/fitc_class_02.py @@ -0,0 +1,66 @@ +import pylab as pb +import numpy as np +import GPy +pb.close('all') + +seed=10000 +"""Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood. + +:param model_type: type of model to fit ['Full', 'FITC', 'DTC']. +:param seed : seed value for data generation. +:type seed: int +:param inducing : number of inducing variables (only used for 'FITC' or 'DTC'). +:type inducing: int +""" + +data = GPy.util.datasets.crescent_data(seed=seed) + +# Kernel object +kernel = GPy.kern.rbf(data['X'].shape[1]) + GPy.kern.white(data['X'].shape[1]) + +# Likelihood object +distribution = GPy.likelihoods.likelihood_functions.probit() +likelihood = GPy.likelihoods.EP(data['Y'],distribution) + +sample = np.random.randint(0,data['X'].shape[0],10) +Z = data['X'][sample,:] +# create sparse GP EP model +#m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z) +m = GPy.models.generalized_FITC(data['X'],likelihood=likelihood,kernel=kernel,Z=Z) +m.ensure_default_constraints() +m.set('len',10.) + +m.update_likelihood_approximation() + +# optimize +m.optimize() +print(m) + +# plot +m.plot() +fitc = m + +pb.figure() +# Kernel object +kernel = GPy.kern.rbf(data['X'].shape[1]) + GPy.kern.white(data['X'].shape[1]) + +# Likelihood object +distribution = GPy.likelihoods.likelihood_functions.probit() +likelihood = GPy.likelihoods.EP(data['Y'],distribution) + +sample = np.random.randint(0,data['X'].shape[0],10) +Z = data['X'][sample,:] +# create sparse GP EP model +m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z) +m.ensure_default_constraints() +m.set('len',10.) + +m.update_likelihood_approximation() + +# optimize +m.optimize() +print(m) + +# plot +m.plot() +variational = m diff --git a/GPy/examples/fitc_regr_01.py b/GPy/examples/fitc_regr_01.py new file mode 100644 index 00000000..815e31a7 --- /dev/null +++ b/GPy/examples/fitc_regr_01.py @@ -0,0 +1,53 @@ +import pylab as pb +import numpy as np +import GPy +pb.ion() +pb.close('all') + +N = 400 +M = 10 +# sample inputs and outputs +X = np.random.uniform(-3.,3.,(N,1)) +Y = np.sin(X)+np.random.randn(N,1)*0.05 + +"""Run a 1D example of a sparse GP regression.""" +""" +rbf = GPy.kern.rbf(1) +noise = GPy.kern.white(1) +kernel = rbf + noise +Z = np.random.uniform(-3,3,(M,1)) +likelihood = GPy.likelihoods.Gaussian(Y) +m = GPy.models.sparse_GP(X, likelihood, kernel, Z) +m.scale_factor=10000 +m.constrain_positive('(variance|lengthscale|precision)') +m.checkgrad(verbose=1) +m.optimize('tnc', messages = 1) +pb.figure() +m.plot() + +variational = m +""" + +# construct kernel +rbf = GPy.kern.rbf(1) +noise = GPy.kern.white(1) +kernel = rbf + noise +#Z = np.random.uniform(-3,3,(M,1)) +Z = variational.Z +likelihood = GPy.likelihoods.Gaussian(Y) +# create simple GP model +m = GPy.models.generalized_FITC(X, likelihood, kernel, Z=Z) +m.constrain_positive('(variance|lengthscale|precision)') +#m.constrain_fixed('iip') +m.checkgrad(verbose=1) +m.optimize('tnc', messages = 1) +#pb.figure() +#m.plot() +""" +Xnew = X.copy().flatten() +Xnew.sort() +Xnew = Xnew[:,None] +mean,var,low,up = m.predict(Xnew) +GPy.util.plot.gpplot(Xnew,mean,low,up) +fitc = m +""" diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py index 739059a0..505f442a 100644 --- a/GPy/models/generalized_FITC.py +++ b/GPy/models/generalized_FITC.py @@ -32,14 +32,8 @@ class generalized_FITC(sparse_GP): def __init__(self, X, likelihood, kernel, Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False): - if type(Z) == int: - #FIXME - self.M = Z - self.Z = (np.random.random_sample(self.D*self.M)*(self.X.max()-self.X.min())+self.X.min()).reshape(self.M,-1) - elif type(Z) == np.ndarray: - self.Z = Z - self.M = self.Z.shape[0] - + self.Z = Z + self.M = self.Z.shape[0] self._precision = likelihood.precision sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False) @@ -64,25 +58,26 @@ class generalized_FITC(sparse_GP): else: self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0) self._precision = self.likelihood.precision # Save the true precision - self.likelihood.precision = self.likelihood.precision/(1. + self.likelihood.precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation + self.likelihood.precision = self._precision/(1. + self._precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation self._set_params(self._get_params()) # update the GP - def _FITC_computations(self): """ FITC approximation doesn't have the correction term in the log-likelihood bound, - but adds a diagonal term to the covariance matrix. + but adds a diagonal term to the covariance matrix: diag(Knn - Qnn). This function: - - computes the diagonal term - - eliminates the extra terms computed in the sparse_GP approximation + - computes the FITC diagonal term + - removes the extra terms computed in the sparse_GP approximation - computes the likelihood gradients wrt the true precision. """ + #NOTE the true precison is now '_precison' not 'precision' if self.likelihood.is_heteroscedastic: + # Compute generalized FITC's diagonal term of the covariance self.Qnn = mdot(self.psi1.T,self.Kmmi,self.psi1) self.Diag0 = self.psi0 - np.diag(self.Qnn) - self.Diag = self.Diag0/(1.+ self.Diag0 * self._precision.flatten()) + self.P = (self.Diag / self.Diag0)[:,None] * self.psi1.T self.RPT0 = np.dot(self.Lmi,self.psi1) self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - self.Diag/(self.Diag0**2))[:,None]*self.RPT0.T)) @@ -92,14 +87,13 @@ class generalized_FITC(sparse_GP): self.w = self.Diag * self.likelihood.v_tilde self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) self.mu = self.w + np.dot(self.P,self.gamma) - self.mu_tilde = (self.likelihood.v_tilde/self.likelihood.tau_tilde)[:,None] # Remove extra term from dL_dpsi1 self.dL_dpsi1 += -mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB else: + raise NotImplementedError, "homoscedastic fitc not implemented" # Remove extra term from dL_dpsi1 - self.dL_dpsi1 += -mdot(self.Kmmi,self.psi1*self.likelihood.precision) #dB - + #self.dL_dpsi1 += -mdot(self.Kmmi,self.psi1*self.likelihood.precision) #dB sf = self.scale_factor sf2 = sf**2 @@ -110,18 +104,16 @@ class generalized_FITC(sparse_GP): #the partial derivative vector for the likelihood if self.likelihood.Nparams == 0: - #save computation here. self.partial_for_likelihood = None elif self.likelihood.is_heteroscedastic: - raise NotImplementedError, "heteroscedatic derivates not implemented" + raise NotImplementedError, "heteroscedastic derivates not implemented" else: + raise NotImplementedError, "homoscedastic derivatives not implemented" #likelihood is not heterscedatic - self.partial_for_likelihood = - 0.5 * self.N*self.D*self._precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self._precision**2 - self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self._precision - self.partial_for_likelihood += self._precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1)) - + #self.partial_for_likelihood = - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self.likelihood.precision**2 + #self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision + #self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1)) #TODO partial derivative vector for the likelihood not implemented - #NOTE the true precison is now '_precison' not 'precision' def dL_dtheta(self): """ @@ -147,23 +139,15 @@ class generalized_FITC(sparse_GP): D = 0.5*np.trace(self.Cpsi1VVpsi1) return A+C+D - def _raw_predict(self, Xnew, slices, full_cov=True): + def _raw_predict(self, Xnew, slices, full_cov=False): if self.likelihood.is_heteroscedastic: """ - Make a prediction for the vsGP model + Make a prediction for the generalized FITC model Arguments --------- X : Input prediction data - Nx1 numpy array (floats) """ - Kx = self.kern.K(self.Z, Xnew) - if full_cov: - Kxx = self.kern.K(Xnew) - else: - Kxx = self.kern.K(Xnew)#FIXME - #raise NotImplementedError - #Kxx = self.kern.Kdiag(Xnew) - # q(u|f) = N(u| R0i*mu_u*f, R0i*C*R0i.T) # Ci = I + (RPT0)Di(RPT0).T @@ -184,13 +168,19 @@ class generalized_FITC(sparse_GP): self.mu_H = mu_H Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T)) # q(f_star|y) = N(f_star|mu_star,sigma2_star) + Kx = self.kern.K(self.Z, Xnew) KR0T = np.dot(Kx.T,self.Lmi.T) mu_star = np.dot(KR0T,mu_H) - sigma2_star = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) - vdiag = np.diag(sigma2_star) - return mu_star[:,None],vdiag[:,None] - else: - """Internal helper function for making predictions, does not account for normalization""" + if full_cov: + Kxx = self.kern.K(Xnew) + var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) + else: + Kxx = self.kern.Kdiag(Xnew) + Kxx_ = self.kern.K(Xnew) + var_ = Kxx_ + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) + var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.M),KR0T.T),0))[:,None] + return mu_star[:,None],var + """ Kx = self.kern.K(self.Z, Xnew) mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) if full_cov: @@ -199,5 +189,19 @@ class generalized_FITC(sparse_GP): else: Kxx = self.kern.Kdiag(Xnew) var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0) - + a = kjk return mu,var[:,None] + """ + else: + raise NotImplementedError, "homoscedastic fitc not implemented" + """ + Kx = self.kern.K(self.Z, Xnew) + mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) + if full_cov: + Kxx = self.kern.K(Xnew) + var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE this won't work for plotting + else: + Kxx = self.kern.Kdiag(Xnew) + var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0) + return mu,var[:,None] + """