From a0df861e8c822f47431e4650437e59a5ebbedb02 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 10 May 2013 17:33:16 +0100 Subject: [PATCH] added swissroll example --- GPy/examples/dimensionality_reduction.py | 65 +++++++++++++++++++++--- 1 file changed, 57 insertions(+), 8 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index e7c775dd..05fb88b5 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -2,13 +2,10 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -import pylab as pb -from matplotlib import pyplot as plt, pyplot +from matplotlib import pyplot as plt import GPy from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM -from GPy.util.datasets import simulation_BGPLVM -from GPy.core.transformations import square, logexp_clipped default_seed = np.random.seed(123344) @@ -47,10 +44,11 @@ def BGPLVM(seed=default_seed): def GPLVM_oil_100(optimize=True): data = GPy.util.datasets.oil_100() + Y = data['X'] # 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 = GPy.models.GPLVM(Y, 6, kernel=kernel) m.data_labels = data['Y'].argmax(axis=1) # optimize @@ -63,22 +61,66 @@ def GPLVM_oil_100(optimize=True): m.plot_latent(labels=m.data_labels) return m -def BGPLVM_oil(optimize=True, N=100, Q=10, M=20, max_f_eval=300, plot=False): +def swiss_roll(optimize=True, N=1000, M=15, Q=4): + from GPy.util.datasets import swiss_roll + from GPy.core.transformations import logexp_clipped + + data = swiss_roll(N=N) + Y = data['Y'] + Y -= Y.mean(0) + Y /= Y.std(0) + + try: + from sklearn.manifold.isomap import Isomap + iso = Isomap().fit(Y) + X = iso.embedding_ + if Q > 2: + X = np.hstack((X, np.random.randn(N, Q - 2))) + except ImportError: + X = np.random.randn(N, Q) + + var = .5 + S = (var * np.ones_like(X) + np.clip(np.random.randn(N, Q) * var ** 2, -(1 - var), (1 - var))) + .001 + Z = np.random.permutation(X)[:M] + + kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, 2) + + m = Bayesian_GPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel) + +# m.constrain('variance|length', logexp_clipped()) + m['lengthscale'] = X.var(0) / X.var(0).max() + m['noise'] = Y.var() / 100. + m.ensure_default_constraints() + + if optimize: + m.optimize('scg', messages=1) + return m + +def BGPLVM_oil(optimize=True, N=100, Q=10, M=10, max_f_eval=1e3, plot=False, **k): data = GPy.util.datasets.oil() + from GPy.core.transformations import logexp_clipped # create simple GP model kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2)) Y = data['X'][:N] - m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M) + Y -= Y.mean(0) + + m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M, **k) m.data_labels = data['Y'][:N].argmax(axis=1) m.constrain('variance', logexp_clipped()) m.constrain('length', logexp_clipped()) - m['lengt'] = 100. + m['lengt'] = 1. + m['noise'] = Y.var() / 100. + m.ensure_default_constraints() # optimize if optimize: + m.unconstrain('X'); m.constrain_fixed('X') + m.optimize('scg', messages=1, max_f_eval=10) + m.unconstrain('X'); m.constrain('X_var', logexp_clipped()) + m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.) m.optimize('scg', messages=1, max_f_eval=150) @@ -115,6 +157,8 @@ def oil_100(): # m.plot_latent(labels=data['Y'].argmax(axis=1)) return m + + def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False): x = np.linspace(0, 4 * np.pi, N)[:, None] s1 = np.vectorize(lambda x: np.sin(x)) @@ -178,6 +222,7 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False): return slist, [S1, S2, S3], Ylist def bgplvm_simulation_matlab_compare(): + from GPy.util.datasets import simulation_BGPLVM sim_data = simulation_BGPLVM() Y = sim_data['Y'] S = sim_data['S'] @@ -213,6 +258,8 @@ def bgplvm_simulation(burnin='scg', plot_sim=False, max_burnin=100, true_X=False, do_opt=True, max_f_eval=1000): + from GPy.core.transformations import logexp_clipped + D1, D2, D3, N, M, Q = 15, 8, 8, 350, 3, 6 slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim) @@ -317,6 +364,8 @@ def mrd_simulation(plot_sim=False): from GPy.models import mrd from GPy import kern + from GPy.core.transformations import logexp_clipped + reload(mrd); reload(kern) # k = kern.rbf(2, ARD=True) + kern.bias(2) + kern.white(2)