From 4953d71b7a3cfe19c29f5155d0285ffb9d7b2025 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 11 Apr 2013 14:54:25 +0100 Subject: [PATCH] first trivial model touches --- GPy/examples/dimensionality_reduction.py | 71 +++++++++++++++++------- GPy/models/GPLVM.py | 4 +- GPy/models/__init__.py | 3 + 3 files changed, 56 insertions(+), 22 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 61a4abd8..a82e0af4 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -6,26 +6,27 @@ import pylab as pb from matplotlib import pyplot as plt import GPy +from GPy.models.mrd import MRD default_seed = np.random.seed(123344) -def BGPLVM(seed = default_seed): +def BGPLVM(seed=default_seed): N = 10 M = 3 Q = 2 D = 4 - #generate GPLVM-like data + # generate GPLVM-like data X = np.random.rand(N, Q) k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) K = k.K(X) - Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y = np.random.multivariate_normal(np.zeros(N), K, D).T - k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q) + k = GPy.kern.linear(Q, ARD=True) + GPy.kern.white(Q) # k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q) # k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) # k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001) - m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) + m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M) m.constrain_positive('(rbf|bias|noise|white|S)') # m.constrain_fixed('S', 1) @@ -38,44 +39,44 @@ def BGPLVM(seed = default_seed): # pb.title('After optimisation') m.ensure_default_constraints() m.randomize() - m.checkgrad(verbose = 1) + m.checkgrad(verbose=1) return m -def GPLVM_oil_100(optimize=True,M=15): +def GPLVM_oil_100(optimize=True, M=15): data = GPy.util.datasets.oil_100() # create simple GP model - kernel = GPy.kern.rbf(6, ARD = True) + GPy.kern.bias(6) + kernel = GPy.kern.rbf(6, ARD=True) + GPy.kern.bias(6) m = GPy.models.GPLVM(data['X'], 6, kernel=kernel, M=M) m.data_labels = data['Y'].argmax(axis=1) # optimize m.ensure_default_constraints() if optimize: - m.optimize('scg',messages=1) + m.optimize('scg', messages=1) # plot print(m) m.plot_latent(labels=m.data_labels) return m -def BGPLVM_oil(optimize=True,N=100,Q=10,M=15): +def BGPLVM_oil(optimize=True, N=100, Q=10, M=15): data = GPy.util.datasets.oil() # create simple GP model - 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) + 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) # optimize if optimize: - m.constrain_fixed('noise',0.05) + m.constrain_fixed('noise', 0.05) m.ensure_default_constraints() - m.optimize('scg',messages=1) + m.optimize('scg', messages=1) m.unconstrain('noise') m.constrain_positive('noise') - m.optimize('scg',messages=1) + m.optimize('scg', messages=1) else: m.ensure_default_constraints() @@ -83,7 +84,7 @@ def BGPLVM_oil(optimize=True,N=100,Q=10,M=15): print(m) m.plot_latent(labels=m.data_labels) pb.figure() - pb.bar(np.arange(m.kern.D),1./m.input_sensitivity()) + pb.bar(np.arange(m.kern.D), 1. / m.input_sensitivity()) return m def oil_100(): @@ -96,7 +97,37 @@ def oil_100(): # plot print(m) - #m.plot_latent(labels=data['Y'].argmax(axis=1)) + # m.plot_latent(labels=data['Y'].argmax(axis=1)) + return m + +def mrd_simulation(): + num = 2 + ard1 = np.array([1., 1, 0, 0], dtype=float) + ard2 = np.array([0., 1, 1, 0], dtype=float) + ard1[ard1 == 0] = 1E+10 + ard2[ard2 == 0] = 1E+10 + + make_params = lambda ard: np.hstack([[1], ard, [1, .3]]) + + D1, D2, N, M, Q = 50, 100, 150, 15, 4 + X = np.random.randn(N, Q) + + k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard1) + GPy.kern.bias(Q, 1) + GPy.kern.white(Q, 0.0001) + Y1 = np.random.multivariate_normal(np.zeros(N), k.K(X), D1).T + Y1 -= Y1.mean(0) + + k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard2) + GPy.kern.bias(Q, 1) + GPy.kern.white(Q, 0.0001) + Y2 = np.random.multivariate_normal(np.zeros(N), k.K(X), D2).T + Y2 -= Y2.mean(0) + + k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q, 1.0) + + m = MRD(Y1, Y2, Q=Q, M=M, kernel=k, _debug=False) + m.ensure_default_constraints() + + m.optimize(messages=1, max_f_eval=5000) + + import ipdb;ipdb.set_trace() return m def brendan_faces(): @@ -109,7 +140,7 @@ def brendan_faces(): m.optimize(messages=1, max_f_eval=10000) ax = m.plot_latent() - y = m.likelihood.Y[0,:] + y = m.likelihood.Y[0, :] data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, invert=False, scale=False) lvm_visualizer = GPy.util.visualize.lvm(m, data_show, ax) raw_input('Press enter to finish') @@ -120,13 +151,13 @@ def brendan_faces(): def stick(): data = GPy.util.datasets.stick() m = GPy.models.GPLVM(data['Y'], 2) - + # optimize m.ensure_default_constraints() m.optimize(messages=1, max_f_eval=10000) ax = m.plot_latent() - y = m.likelihood.Y[0,:] + y = m.likelihood.Y[0, :] data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect']) lvm_visualizer = GPy.util.visualize.lvm(m, data_show, ax) raw_input('Press enter to finish') diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py index cc4be70e..2b1002eb 100644 --- a/GPy/models/GPLVM.py +++ b/GPy/models/GPLVM.py @@ -55,7 +55,7 @@ class GPLVM(GP): def plot(self): assert self.likelihood.Y.shape[1]==2 - pb.scatter(self.likelihood.Y[:,0],self.likelihood.Y[:,1],40,self.X[:,0].copy(),linewidth=0,cmap=pb.cm.jet) + pb.scatter(self.likelihood.Y[:,0],self.likelihood.Y[:,1],40,self.X[:,0].copy(),linewidth=0,cmap=pb.cm.jet) # @UndefinedVariable Xnew = np.linspace(self.X.min(),self.X.max(),200)[:,None] mu, var, upper, lower = self.predict(Xnew) pb.plot(mu[:,0], mu[:,1],'k',linewidth=1.5) @@ -90,7 +90,7 @@ class GPLVM(GP): 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') + pb.imshow(var.reshape(resolution,resolution).T[::-1,:],extent=[xmin[0],xmax[0],xmin[1],xmax[1]],cmap=pb.cm.binary,interpolation='bilinear') # @UndefinedVariable for i,ul in enumerate(np.unique(labels)): diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py index f442dc67..91cc60e3 100644 --- a/GPy/models/__init__.py +++ b/GPy/models/__init__.py @@ -11,4 +11,7 @@ from warped_GP import warpedGP from sparse_GPLVM import sparse_GPLVM from uncollapsed_sparse_GP import uncollapsed_sparse_GP from Bayesian_GPLVM import Bayesian_GPLVM +import mrd +MRD = mrd.MRD +del mrd from generalized_FITC import generalized_FITC