diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index ac39fd66..ef6e37fb 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -7,26 +7,27 @@ from matplotlib import pyplot as plt, cm import GPy from GPy.core.transformations import logexp from GPy.models.bayesian_gplvm import BayesianGPLVM +from GPy.likelihoods.gaussian import Gaussian default_seed = np.random.seed(123344) def BGPLVM(seed=default_seed): N = 10 num_inducing = 3 - Q = 2 - D = 4 + Q = 5 + D = 10 # generate GPLVM-like data X = np.random.rand(N, Q) - k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.01) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N), K, Q).T + lik = Gaussian(Y, normalize=True) - k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) + GPy.kern.white(Q) - # k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q) + k = GPy.kern.rbf_inv(Q) + GPy.kern.bias(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.BayesianGPLVM(Y, Q, kernel=k, num_inducing=num_inducing) + m = GPy.models.BayesianGPLVM(lik, Q, kernel=k, num_inducing=num_inducing) # m.constrain_positive('(rbf|bias|noise|white|S)') # m.constrain_fixed('S', 1) @@ -70,16 +71,16 @@ def sparseGPLVM_oil(optimize=True, N=100, Q=6, num_inducing=15, max_iters=50): # create simple GP model kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) - m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing = num_inducing) + m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing) m.data_labels = data['Y'].argmax(axis=1) # optimize if optimize: - m.optimize('scg', messages=1, max_iters = max_iters) + m.optimize('scg', messages=1, max_iters=max_iters) # plot print(m) - #m.plot_latent(labels=m.data_labels) + # m.plot_latent(labels=m.data_labels) return m def swiss_roll(optimize=True, N=1000, num_inducing=15, Q=4, sigma=.2, plot=False): @@ -136,12 +137,12 @@ def swiss_roll(optimize=True, N=1000, num_inducing=15, Q=4, sigma=.2, plot=False m.optimize('scg', messages=1) return m -def BGPLVM_oil(optimize=True, N=200, Q=10, num_inducing=15, max_iters=50, plot=False, **k): +def BGPLVM_oil(optimize=True, N=200, Q=10, num_inducing=15, max_iters=150, plot=False, **k): np.random.seed(0) data = GPy.util.datasets.oil() # 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)) + kernel = GPy.kern.rbf_inv(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2)) Y = data['X'][:N] Yn = Y - Y.mean(0) Yn /= Yn.std(0) @@ -150,14 +151,14 @@ def BGPLVM_oil(optimize=True, N=200, Q=10, num_inducing=15, max_iters=50, plot=F m.data_labels = data['Y'][:N].argmax(axis=1) # m.constrain('variance|leng', logexp_clipped()) - m['.*lengt'] = 1. # m.X.var(0).max() / m.X.var(0) + m['.*lengt'] = m.X.var(0).max() / m.X.var(0) m['noise'] = Yn.var() / 100. # optimize if optimize: m.constrain_fixed('noise') - m.optimize('scg', messages=1, max_iters=100, gtol=.05) + m.optimize('scg', messages=1, max_iters=200, gtol=.05) m.constrain_positive('noise') m.optimize('scg', messages=1, max_iters=max_iters, gtol=.05) @@ -208,7 +209,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): Y1 += .3 * np.random.randn(*Y1.shape) Y2 += .2 * np.random.randn(*Y2.shape) - Y3 += .1 * np.random.randn(*Y3.shape) + Y3 += .25 * np.random.randn(*Y3.shape) Y1 -= Y1.mean(0) Y2 -= Y2.mean(0) @@ -263,16 +264,16 @@ def bgplvm_simulation_matlab_compare(): def bgplvm_simulation(optimize='scg', plot=True, - max_iters=2e4): + max_iters=2e4, + plot_sim=False): # from GPy.core.transformations import logexp_clipped - D1, D2, D3, N, num_inducing, Q = 15, 8, 8, 100, 3, 5 - slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot) + D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 300, 23, 6 + slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) from GPy.models import mrd from GPy import kern reload(mrd); reload(kern) - Y = Ylist[0] k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q) @@ -280,7 +281,6 @@ def bgplvm_simulation(optimize='scg', # m.constrain('variance|noise', logexp_clipped()) m['noise'] = Y.var() / 100. - m['linear_variance'] = .01 if optimize: print "Optimizing model:" @@ -346,7 +346,7 @@ def brendan_faces(): def stick_play(range=None, frame_rate=15): data = GPy.util.datasets.stick() # optimize - if range==None: + if range == None: Y = data['Y'].copy() else: Y = data['Y'][range[0]:range[1], :].copy() @@ -374,10 +374,10 @@ def stick_bgplvm(model=None): data = GPy.util.datasets.stick() Q = 6 kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2)) - m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20,kernel=kernel) + m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel) # optimize m.ensure_default_constraints() - m.optimize(messages=1, max_f_eval=3000,xtol=1e-300,ftol=1e-300) + m.optimize(messages=1, max_f_eval=3000, xtol=1e-300, ftol=1e-300) m._set_params(m._get_params()) plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2) plt.sca(latent_axes)