diff --git a/GPy/core/model.py b/GPy/core/model.py index 9216aea6..645f7228 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -423,6 +423,28 @@ class model(parameterised): grad_string = "{0:^{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name,r,d,g, ng, c0 = cols[0]+9, c1 = cols[1], c2 = cols[2], c3 = cols[3], c4 = cols[4]) print grad_string + def input_sensitivity(self): + """ + return an array describing the sesitivity of the model to each input + + NB. Right now, we're basing this on the lengthscales (or variances) of the kernel. + TODO: proper sensitivity analysis + """ + + if not hasattr(self,'kern'): + raise ValueError, "this model has no kernel" + + k = [p for p in self.kern.parts if p.name in ['rbf','linear']] + if (not len(k)==1) or (not k[0].ARD): + raise ValueError, "cannot determine sensitivity for this kernel" + k = k[0] + + if k.name=='rbf': + return k.lengthscale + elif k.name=='linear': + return 1./k.variances + + def pseudo_EM(self,epsilon=.1,**kwargs): """ TODO: Should this not bein the GP class? @@ -469,3 +491,4 @@ class model(parameterised): iteration += 1 if stop: print "%s iterations." %iteration + diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index fed79417..b3509ec5 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -47,7 +47,7 @@ def GPLVM_oil_100(optimize=True,M=15): # create simple GP model kernel = GPy.kern.rbf(6, ARD = True) + GPy.kern.bias(6) - m = GPy.models.Bayesian_GPLVM(data['X'], 6, kernel=kernel, M=M) + m = GPy.models.GPLVM(data['X'], 6, kernel=kernel, M=M) m.data_labels = data['Y'].argmax(axis=1) # optimize @@ -60,22 +60,32 @@ def GPLVM_oil_100(optimize=True,M=15): m.plot_latent(labels=m.data_labels) return m -def BGPLVM_oil_100(optimize=True): - data = GPy.util.datasets.oil_100() +def BGPLVM_oil(optimize=True,N=100,Q=10,M=15): + data = GPy.util.datasets.oil() # create simple GP model - kernel = GPy.kern.rbf(6, ARD = True) + GPy.kern.bias(6) - m = GPy.models.GPLVM(data['X'], 6, kernel = kernel) - m.data_labels = data['Y'].argmax(axis=1) + kernel = GPy.kern.rbf(Q, ARD = True) + GPy.kern.bias(Q) + GPy.kern.white(Q,0.001) + m = GPy.models.Bayesian_GPLVM(data['X'][:N], Q, kernel = kernel,M=M) + m.data_labels = data['Y'][:N].argmax(axis=1) + + #initial conditions # optimize - m.ensure_default_constraints() if optimize: - m.optimize(messages=1) + m.constrain_fixed('noise',0.05) + m.ensure_default_constraints() + m.optimize('scg',messages=1) + m.unconstrain('noise') + m.constrain_positive('noise') + m.optimize('scg',messages=1) + else: + m.ensure_default_constraints() # plot print(m) - m.plot_latent(labels=data['Y'].argmax(axis=1)) + m.plot_latent(labels=m.data_labels) + pb.figure() + pb.bar(np.arange(m.kern.D),1./m.input_sensitivity()) return m def oil_100(): diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index 8f9759c3..0f5758ba 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -22,12 +22,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): :type init: 'PCA'|'random' """ - def __init__(self, Y, Q, X = None, S = None, init='PCA', M=10, Z=None, kernel=None, **kwargs): + def __init__(self, Y, Q, X = None, X_uncertainty = None, init='PCA', M=10, Z=None, kernel=None, **kwargs): if X == None: X = self.initialise_latent(init, Q, Y) - if S is None: - S = np.ones_like(X) * 1e-2# + if X_uncertainty is None: + X_uncertainty = np.ones_like(X) * 0.5 if Z is None: Z = np.random.permutation(X.copy())[:M] @@ -37,7 +37,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): kernel = kern.rbf(Q) + kern.white(Q) - sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_uncertainty=S, **kwargs) + sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_uncertainty=X_uncertainty, **kwargs) def _get_param_names(self): X_names = sum([['X_%i_%i'%(n,q) for q in range(self.Q)] for n in range(self.N)],[]) @@ -84,6 +84,14 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): def _log_likelihood_gradients(self): return np.hstack((self.dL_dmuS().flatten(), sparse_GP._log_likelihood_gradients(self))) - def plot_latent(self, *args, **kwargs): - input_1, input_2 = GPLVM.plot_latent(self, *args, **kwargs) + def plot_latent(self, which_indices=None,*args, **kwargs): + + if which_indices is None: + try: + input_1, input_2 = np.argsort(self.input_sensitivity())[:2] + except: + raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'" + else: + input_1, input_2 = which_indices + GPLVM.plot_latent(self, which_indices=[input_1, input_2],*args, **kwargs) pb.plot(self.Z[:, input_1], self.Z[:, input_2], '^w') diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py index 03e9b715..cc4be70e 100644 --- a/GPy/models/GPLVM.py +++ b/GPy/models/GPLVM.py @@ -67,7 +67,7 @@ class GPLVM(GP): """ util.plot.Tango.reset() - + if labels is None: labels = np.ones(self.N) if which_indices is None: @@ -77,23 +77,19 @@ class GPLVM(GP): if self.Q==2: input_1, input_2 = 0,1 else: - #try to find a linear of RBF kern in the kernel - k = [p for p in self.kern.parts if p.name in ['rbf','linear']] - if (not len(k)==1) or (not k[0].ARD): + try: + input_1, input_2 = np.argsort(self.input_sensitivity())[:2] + except: raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'" - input_1, input_2 = self.lengthscale_order() - k = k[0] - if k.name=='rbf': - input_1, input_2 = np.argsort(k.lengthscale)[:2] - elif k.name=='linear': - input_1, input_2 = np.argsort(k.variances)[::-1][:2] + else: + input_1, input_2 = which_indices #first, plot the output variance as a function of the latent space Xtest, xx,yy,xmin,xmax = util.plot.x_frame2D(self.X[:,[input_1, input_2]],resolution=resolution) - Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1])) - Xtest_full[:, :2] = Xtest - mu, var, low, up = self.predict(Xtest_full) - var = var.mean(axis=1) # this was var[:, :2] edit by Neil + Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1])) + Xtest_full[:, :2] = Xtest + mu, var, low, up = self.predict(Xtest_full) + var = var.mean(axis=1) # this was var[:, :2] edit by Neil pb.imshow(var.reshape(resolution,resolution).T[::-1,:],extent=[xmin[0],xmax[0],xmin[1],xmax[1]],cmap=pb.cm.binary,interpolation='bilinear')