diff --git a/GPy/examples/classification.py b/GPy/examples/classification.py index af1c8c7a..ce3bd6d8 100644 --- a/GPy/examples/classification.py +++ b/GPy/examples/classification.py @@ -7,39 +7,47 @@ import GPy default_seed = 10000 + def oil(num_inducing=50, max_iters=100, kernel=None, optimize=True, plot=True): """ Run a Gaussian process classification on the three phase oil data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood. """ - try:import pods - except ImportError:raise ImportWarning('Need pods for example datasets. See https://github.com/sods/ods, or pip install pods.') + try: + import pods + except ImportError: + raise ImportWarning( + "Need pods for example datasets. See https://github.com/sods/ods, or pip install pods." + ) data = pods.datasets.oil() - X = data['X'] - Xtest = data['Xtest'] - Y = data['Y'][:, 0:1] - Ytest = data['Ytest'][:, 0:1] - Y[Y.flatten()==-1] = 0 - Ytest[Ytest.flatten()==-1] = 0 + X = data["X"] + Xtest = data["Xtest"] + Y = data["Y"][:, 0:1] + Ytest = data["Ytest"][:, 0:1] + Y[Y.flatten() == -1] = 0 + Ytest[Ytest.flatten() == -1] = 0 # Create GP model - m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, num_inducing=num_inducing) + m = GPy.models.SparseGPClassification( + X, Y, kernel=kernel, num_inducing=num_inducing + ) m.Ytest = Ytest # Contrain all parameters to be positive - #m.tie_params('.*len') - m['.*len'] = 10. + # m.tie_params('.*len') + m[".*len"] = 10.0 # Optimize if optimize: m.optimize(messages=1) print(m) - #Test + # Test probs = m.predict(Xtest)[0] GPy.util.classification.conf_matrix(probs, Ytest) return m + def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True): """ Simple 1D classification example using EP approximation @@ -48,26 +56,31 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True): :type seed: int """ - try:import pods - except ImportError:raise ImportWarning('Need pods for example datasets. See https://github.com/sods/ods, or pip install pods.') + try: + import pods + except ImportError: + raise ImportWarning( + "Need pods for example datasets. See https://github.com/sods/ods, or pip install pods." + ) data = pods.datasets.toy_linear_1d_classification(seed=seed) - Y = data['Y'][:, 0:1] + Y = data["Y"][:, 0:1] Y[Y.flatten() == -1] = 0 # Model definition - m = GPy.models.GPClassification(data['X'], Y) + m = GPy.models.GPClassification(data["X"], Y) # Optimize if optimize: - #m.update_likelihood_approximation() + # m.update_likelihood_approximation() # Parameters optimization: m.optimize() - #m.update_likelihood_approximation() - #m.pseudo_EM() + # m.update_likelihood_approximation() + # m.pseudo_EM() # Plot if plot: from matplotlib import pyplot as plt + fig, axes = plt.subplots(2, 1) m.plot_f(ax=axes[0]) m.plot(ax=axes[1]) @@ -75,6 +88,7 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True): print(m) return m + def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=True): """ Simple 1D classification example using Laplace approximation @@ -84,10 +98,12 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot= """ - try:import pods - except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets') + try: + import pods + except ImportError: + print("pods unavailable, see https://github.com/sods/ods for example datasets") data = pods.datasets.toy_linear_1d_classification(seed=seed) - Y = data['Y'][:, 0:1] + Y = data["Y"][:, 0:1] Y[Y.flatten() == -1] = 0 likelihood = GPy.likelihoods.Bernoulli() @@ -95,18 +111,21 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot= kernel = GPy.kern.RBF(1) # Model definition - m = GPy.core.GP(data['X'], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf) + m = GPy.core.GP( + data["X"], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf + ) # Optimize if optimize: try: - m.optimize('scg', messages=1) + m.optimize("scg", messages=1) except Exception as e: return m # Plot if plot: from matplotlib import pyplot as plt + fig, axes = plt.subplots(2, 1) m.plot_f(ax=axes[0]) m.plot(ax=axes[1]) @@ -114,7 +133,10 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot= print(m) return m -def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, optimize=True, plot=True): + +def sparse_toy_linear_1d_classification( + num_inducing=10, seed=default_seed, optimize=True, plot=True +): """ Sparse 1D classification example @@ -123,15 +145,17 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti """ - try:import pods - except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets') + try: + import pods + except ImportError: + print("pods unavailable, see https://github.com/sods/ods for example datasets") data = pods.datasets.toy_linear_1d_classification(seed=seed) - Y = data['Y'][:, 0:1] + Y = data["Y"][:, 0:1] Y[Y.flatten() == -1] = 0 # Model definition - m = GPy.models.SparseGPClassification(data['X'], Y, num_inducing=num_inducing) - m['.*len'] = 4. + m = GPy.models.SparseGPClassification(data["X"], Y, num_inducing=num_inducing) + m[".*len"] = 4.0 # Optimize if optimize: @@ -140,6 +164,7 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti # Plot if plot: from matplotlib import pyplot as plt + fig, axes = plt.subplots(2, 1) m.plot_f(ax=axes[0]) m.plot(ax=axes[1]) @@ -147,7 +172,10 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti print(m) return m -def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=default_seed, optimize=True, plot=True): + +def sparse_toy_linear_1d_classification_uncertain_input( + num_inducing=10, seed=default_seed, optimize=True, plot=True +): """ Sparse 1D classification example @@ -156,18 +184,23 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de """ - try:import pods - except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets') + try: + import pods + except ImportError: + print("pods unavailable, see https://github.com/sods/ods for example datasets") import numpy as np + data = pods.datasets.toy_linear_1d_classification(seed=seed) - Y = data['Y'][:, 0:1] + Y = data["Y"][:, 0:1] Y[Y.flatten() == -1] = 0 - X = data['X'] - X_var = np.random.uniform(0.3,0.5,X.shape) + X = data["X"] + X_var = np.random.uniform(0.3, 0.5, X.shape) # Model definition - m = GPy.models.SparseGPClassificationUncertainInput(X, X_var, Y, num_inducing=num_inducing) - m['.*len'] = 4. + m = GPy.models.SparseGPClassificationUncertainInput( + X, X_var, Y, num_inducing=num_inducing + ) + m[".*len"] = 4.0 # Optimize if optimize: @@ -176,6 +209,7 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de # Plot if plot: from matplotlib import pyplot as plt + fig, axes = plt.subplots(2, 1) m.plot_f(ax=axes[0]) m.plot(ax=axes[1]) @@ -183,6 +217,7 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de print(m) return m + def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True): """ Simple 1D classification example using a heavy side gp transformation @@ -192,29 +227,41 @@ def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True): """ - try:import pods - except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets') + try: + import pods + except ImportError: + print("pods unavailable, see https://github.com/sods/ods for example datasets") data = pods.datasets.toy_linear_1d_classification(seed=seed) - Y = data['Y'][:, 0:1] + Y = data["Y"][:, 0:1] Y[Y.flatten() == -1] = 0 # Model definition kernel = GPy.kern.RBF(1) - likelihood = GPy.likelihoods.Bernoulli(gp_link=GPy.likelihoods.link_functions.Heaviside()) + likelihood = GPy.likelihoods.Bernoulli( + gp_link=GPy.likelihoods.link_functions.Heaviside() + ) ep = GPy.inference.latent_function_inference.expectation_propagation.EP() - m = GPy.core.GP(X=data['X'], Y=Y, kernel=kernel, likelihood=likelihood, inference_method=ep, name='gp_classification_heaviside') - #m = GPy.models.GPClassification(data['X'], likelihood=likelihood) + m = GPy.core.GP( + X=data["X"], + Y=Y, + kernel=kernel, + likelihood=likelihood, + inference_method=ep, + name="gp_classification_heaviside", + ) + # m = GPy.models.GPClassification(data['X'], likelihood=likelihood) # Optimize if optimize: # Parameters optimization: for _ in range(5): - m.optimize(max_iters=int(max_iters/5)) + m.optimize(max_iters=int(max_iters / 5)) print(m) # Plot if plot: from matplotlib import pyplot as plt + fig, axes = plt.subplots(2, 1) m.plot_f(ax=axes[0]) m.plot(ax=axes[1]) @@ -222,7 +269,15 @@ def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True): print(m) return m -def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=None, optimize=True, plot=True): + +def crescent_data( + model_type="Full", + num_inducing=10, + seed=default_seed, + kernel=None, + optimize=True, + plot=True, +): """ Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood. @@ -234,22 +289,28 @@ def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel= :param kernel: kernel to use in the model :type kernel: a GPy kernel """ - try:import pods - except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets') + try: + import pods + except ImportError: + print("pods unavailable, see https://github.com/sods/ods for example datasets") data = pods.datasets.crescent_data(seed=seed) - Y = data['Y'] - Y[Y.flatten()==-1] = 0 + Y = data["Y"] + Y[Y.flatten() == -1] = 0 - if model_type == 'Full': - m = GPy.models.GPClassification(data['X'], Y, kernel=kernel) + if model_type == "Full": + m = GPy.models.GPClassification(data["X"], Y, kernel=kernel) - elif model_type == 'DTC': - m = GPy.models.SparseGPClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing) - m['.*len'] = 10. + elif model_type == "DTC": + m = GPy.models.SparseGPClassification( + data["X"], Y, kernel=kernel, num_inducing=num_inducing + ) + m[".*len"] = 10.0 - elif model_type == 'FITC': - m = GPy.models.FITCClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing) - m['.*len'] = 3. + elif model_type == "FITC": + m = GPy.models.FITCClassification( + data["X"], Y, kernel=kernel, num_inducing=num_inducing + ) + m[".*len"] = 3.0 if optimize: m.optimize(messages=1) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 9218a0e8..74df5235 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -1,10 +1,12 @@ # Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as _np + default_seed = 123344 # default_seed = _np.random.seed(123344) + def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan=False): """ model for testing purposes. Samples from a GP with rbf kernel and learns @@ -24,7 +26,7 @@ def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan # generate GPLVM-like data X = _np.random.rand(num_inputs, input_dim) lengthscales = _np.random.rand(input_dim) - k = GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True) + k = GPy.kern.RBF(input_dim, 0.5, lengthscales, ARD=True) K = k.K(X) Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, (output_dim,)).T @@ -35,88 +37,102 @@ def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan # k = GPy.kern.RBF(input_dim, .5, 2., ARD=0) + GPy.kern.RBF(input_dim, .3, .2, ARD=0) # k = GPy.kern.RBF(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.linear(input_dim, _np.ones(input_dim) * .2, ARD=True) - p = .3 + p = 0.3 m = GPy.models.BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing) if nan: - m.inference_method = GPy.inference.latent_function_inference.var_dtc.VarDTCMissingData() + m.inference_method = ( + GPy.inference.latent_function_inference.var_dtc.VarDTCMissingData() + ) m.Y[_np.random.binomial(1, p, size=(Y.shape)).astype(bool)] = _np.nan m.parameters_changed() - #=========================================================================== + # =========================================================================== # randomly obstruct data with percentage p - #=========================================================================== + # =========================================================================== # m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing) # m.lengthscales = lengthscales if plot: import matplotlib.pyplot as pb + m.plot() - pb.title('PCA initialisation') + pb.title("PCA initialisation") # m2.plot() # pb.title('PCA initialisation') if optimize: - m.optimize('scg', messages=verbose) + m.optimize("scg", messages=verbose) # m2.optimize('scg', messages=verbose) if plot: m.plot() - pb.title('After optimisation') + pb.title("After optimisation") # m2.plot() # pb.title('After optimisation') return m + def gplvm_oil_100(optimize=True, verbose=1, plot=True): import GPy import pods + data = pods.datasets.oil_100() - Y = data['X'] + Y = data["X"] # create simple GP model kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6) m = GPy.models.GPLVM(Y, 6, kernel=kernel) - m.data_labels = data['Y'].argmax(axis=1) - if optimize: m.optimize('scg', messages=verbose) + m.data_labels = data["Y"].argmax(axis=1) + if optimize: + m.optimize("scg", messages=verbose) if plot: m.plot_latent(labels=m.data_labels) return m -def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50): + +def sparse_gplvm_oil( + optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50 +): import GPy import pods _np.random.seed(0) data = pods.datasets.oil() - Y = data['X'][:N] + Y = data["X"][:N] Y = Y - Y.mean(0) Y /= Y.std(0) # Create the 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.data_labels = data['Y'][:N].argmax(axis=1) + m.data_labels = data["Y"][:N].argmax(axis=1) - if optimize: m.optimize('scg', messages=verbose, max_iters=max_iters) + if optimize: + m.optimize("scg", messages=verbose, max_iters=max_iters) if plot: m.plot_latent(labels=m.data_labels) m.kern.plot_ARD() return m -def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2): + +def swiss_roll( + optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=0.2 +): import GPy from pods.datasets import swiss_roll_generated from GPy.models import BayesianGPLVM data = swiss_roll_generated(num_samples=N, sigma=sigma) - Y = data['Y'] + Y = data["Y"] Y -= Y.mean() Y /= Y.std() - t = data['t'] - c = data['colors'] + t = data["t"] + c = data["colors"] try: from sklearn.manifold.isomap import Isomap + iso = Isomap().fit(Y) X = iso.embedding_ if Q > 2: @@ -127,8 +143,9 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4 if plot: import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # @UnusedImport + fig = plt.figure("Swiss Roll Data") - ax = fig.add_subplot(121, projection='3d') + ax = fig.add_subplot(121, projection="3d") ax.scatter(*Y.T, c=c) ax.set_title("Swiss Roll") @@ -136,60 +153,96 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4 ax.scatter(*X.T[:2], c=c) ax.set_title("BGPLVM init") - var = .5 - S = (var * _np.ones_like(X) + _np.clip(_np.random.randn(N, Q) * var ** 2, - - (1 - var), - (1 - var))) + .001 + var = 0.5 + S = ( + var * _np.ones_like(X) + + _np.clip(_np.random.randn(N, Q) * var ** 2, -(1 - var), (1 - var)) + ) + 0.001 Z = _np.random.permutation(X)[:num_inducing] - 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(Q, ARD=True) + + GPy.kern.Bias(Q, _np.exp(-2)) + + GPy.kern.White(Q, _np.exp(-2)) + ) - m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel) + m = BayesianGPLVM( + Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel + ) m.data_colors = c m.data_t = t if optimize: - m.optimize('bfgs', messages=verbose, max_iters=2e3) + m.optimize("bfgs", messages=verbose, max_iters=2e3) if plot: - fig = plt.figure('fitted') + fig = plt.figure("fitted") ax = fig.add_subplot(111) s = m.input_sensitivity().argsort()[::-1][:2] ax.scatter(*m.X.mean.T[s], c=c) return m -def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k): + +def bgplvm_oil( + optimize=True, + verbose=1, + plot=True, + N=200, + Q=7, + num_inducing=40, + max_iters=1000, + **k +): import GPy from matplotlib import pyplot as plt import numpy as np + _np.random.seed(0) try: import pods + data = pods.datasets.oil() except ImportError: data = GPy.util.datasets.oil() - - kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2)) - Y = data['X'][:N] + kernel = GPy.kern.RBF( + Q, 1.0, 1.0 / _np.random.uniform(0, 1, (Q,)), ARD=True + ) # + GPy.kern.Bias(Q, _np.exp(-2)) + Y = data["X"][:N] m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) - m.data_labels = data['Y'][:N].argmax(axis=1) + m.data_labels = data["Y"][:N].argmax(axis=1) if optimize: - m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: fig, (latent_axes, sense_axes) = plt.subplots(1, 2) m.plot_latent(ax=latent_axes, labels=m.data_labels) data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) - lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable - m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels) - input('Press enter to finish') + lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect( + m.X.mean.values[0:1, :], # @UnusedVariable + m, + data_show, + latent_axes=latent_axes, + sense_axes=sense_axes, + labels=m.data_labels, + ) + input("Press enter to finish") plt.close(fig) return m -def ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k): + +def ssgplvm_oil( + optimize=True, + verbose=1, + plot=True, + N=200, + Q=7, + num_inducing=40, + max_iters=1000, + **k +): import GPy from matplotlib import pyplot as plt import pods @@ -197,39 +250,57 @@ def ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40 _np.random.seed(0) data = pods.datasets.oil() - kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2)) - Y = data['X'][:N] + kernel = GPy.kern.RBF( + Q, 1.0, 1.0 / _np.random.uniform(0, 1, (Q,)), ARD=True + ) # + GPy.kern.Bias(Q, _np.exp(-2)) + Y = data["X"][:N] m = GPy.models.SSGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) - m.data_labels = data['Y'][:N].argmax(axis=1) + m.data_labels = data["Y"][:N].argmax(axis=1) if optimize: - m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: fig, (latent_axes, sense_axes) = plt.subplots(1, 2) m.plot_latent(ax=latent_axes, labels=m.data_labels) data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) - lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable - m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels) - input('Press enter to finish') + lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect( + m.X.mean.values[0:1, :], # @UnusedVariable + m, + data_show, + latent_axes=latent_axes, + sense_axes=sense_axes, + labels=m.data_labels, + ) + input("Press enter to finish") plt.close(fig) return m + def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False): """Simulate some data drawn from a matern covariance and a periodic exponential for use in MRD demos.""" Q_signal = 4 import GPy import numpy as np + np.random.seed(3000) - k = GPy.kern.Matern32(Q_signal, 1., lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1) + k = GPy.kern.Matern32( + Q_signal, 1.0, lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1 + ) for i in range(Q_signal): - k += GPy.kern.PeriodicExponential(1, variance=1., active_dims=[i], period=3., lower=-2, upper=6) + k += GPy.kern.PeriodicExponential( + 1, variance=1.0, active_dims=[i], period=3.0, lower=-2, upper=6 + ) t = np.c_[[np.linspace(-1, 5, N) for _ in range(Q_signal)]].T K = k.K(t) - s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:, :, None] + s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[ + :, :, None + ] - Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS) + Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output( + D1, D2, D3, s1, s2, s3, sS + ) slist = [sS, s1, s2, s3] slist_names = ["sS", "s1", "s2", "s3"] @@ -239,6 +310,7 @@ def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False): from matplotlib import pyplot as plt import matplotlib.cm as cm import itertools + fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) fig.clf() ax = fig.add_subplot(2, 1, 1) @@ -248,13 +320,14 @@ def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False): ax.legend() for i, Y in enumerate(Ylist): ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) - ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable + ax.imshow(Y, aspect="auto", cmap=cm.gray) # @UndefinedVariable ax.set_title("Y{}".format(i + 1)) plt.draw() plt.tight_layout() return slist, [S1, S2, S3], Ylist + def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): """Simulate some data drawn from sine and cosine for use in demos of MRD""" _np.random.seed(1234) @@ -262,7 +335,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): x = _np.linspace(0, 4 * _np.pi, N)[:, None] s1 = _np.vectorize(lambda x: _np.sin(x)) s2 = _np.vectorize(lambda x: _np.cos(x)) - s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) + s3 = _np.vectorize(lambda x: -_np.exp(-_np.cos(2 * x))) sS = _np.vectorize(lambda x: _np.cos(x)) s1 = s1(x) @@ -270,12 +343,18 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): s3 = s3(x) sS = sS(x) - s1 -= s1.mean(); s1 /= s1.std(0) - s2 -= s2.mean(); s2 /= s2.std(0) - s3 -= s3.mean(); s3 /= s3.std(0) - sS -= sS.mean(); sS /= sS.std(0) + s1 -= s1.mean() + s1 /= s1.std(0) + s2 -= s2.mean() + s2 /= s2.std(0) + s3 -= s3.mean() + s3 /= s3.std(0) + sS -= sS.mean() + sS /= sS.std(0) - Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS) + Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output( + D1, D2, D3, s1, s2, s3, sS + ) slist = [sS, s1, s2, s3] slist_names = ["sS", "s1", "s2", "s3"] @@ -285,6 +364,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): from matplotlib import pyplot as plt import matplotlib.cm as cm import itertools + fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) fig.clf() ax = fig.add_subplot(2, 1, 1) @@ -294,13 +374,14 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): ax.legend() for i, Y in enumerate(Ylist): ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) - ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable + ax.imshow(Y, aspect="auto", cmap=cm.gray) # @UndefinedVariable ax.set_title("Y{}".format(i + 1)) plt.draw() plt.tight_layout() return slist, [S1, S2, S3], Ylist + def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS): S1 = _np.hstack([s1, sS]) S2 = _np.hstack([sS]) @@ -308,9 +389,9 @@ def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS): Y1 = S1.dot(_np.random.randn(S1.shape[1], D1)) Y2 = S2.dot(_np.random.randn(S2.shape[1], D2)) Y3 = S3.dot(_np.random.randn(S3.shape[1], D3)) - Y1 += .3 * _np.random.randn(*Y1.shape) - Y2 += .2 * _np.random.randn(*Y2.shape) - Y3 += .25 * _np.random.randn(*Y3.shape) + Y1 += 0.3 * _np.random.randn(*Y1.shape) + Y2 += 0.2 * _np.random.randn(*Y2.shape) + Y3 += 0.25 * _np.random.randn(*Y3.shape) Y1 -= Y1.mean(0) Y2 -= Y2.mean(0) Y3 -= Y3.mean(0) @@ -319,10 +400,10 @@ def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS): Y3 /= Y3.std(0) return Y1, Y2, Y3, S1, S2, S3 -def bgplvm_simulation(optimize=True, verbose=1, - plot=True, plot_sim=False, - max_iters=2e4, - ): + +def bgplvm_simulation( + optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, +): from GPy import kern from GPy.models import BayesianGPLVM @@ -332,22 +413,21 @@ def bgplvm_simulation(optimize=True, verbose=1, k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k) - m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape) - m.likelihood.variance = .1 + m.X.variance[:] = _np.random.uniform(0, 0.01, m.X.shape) + m.likelihood.variance = 0.1 if optimize: print("Optimizing model:") - m.optimize('bfgs', messages=verbose, max_iters=max_iters, - gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m -def gplvm_simulation(optimize=True, verbose=1, - plot=True, plot_sim=False, - max_iters=2e4, - ): + +def gplvm_simulation( + optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, +): from GPy import kern from GPy.models import GPLVM @@ -357,20 +437,20 @@ def gplvm_simulation(optimize=True, verbose=1, k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) m = GPLVM(Y, Q, init="PCA", kernel=k) - m.likelihood.variance = .1 + m.likelihood.variance = 0.1 if optimize: print("Optimizing model:") - m.optimize('bfgs', messages=verbose, max_iters=max_iters, - gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m -def ssgplvm_simulation(optimize=True, verbose=1, - plot=True, plot_sim=False, - max_iters=2e4, useGPU=False - ): + + +def ssgplvm_simulation( + optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, useGPU=False +): from GPy import kern from GPy.models import SSGPLVM @@ -379,23 +459,30 @@ def ssgplvm_simulation(optimize=True, verbose=1, Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) - m = SSGPLVM(Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True) - m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape) - m.likelihood.variance = .01 + m = SSGPLVM( + Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True + ) + m.X.variance[:] = _np.random.uniform(0, 0.01, m.X.shape) + m.likelihood.variance = 0.01 if optimize: print("Optimizing model:") - m.optimize('bfgs', messages=verbose, max_iters=max_iters, - gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("SSGPLVM Latent Space 1D") m.kern.plot_ARD() return m -def bgplvm_simulation_missing_data(optimize=True, verbose=1, - plot=True, plot_sim=False, - max_iters=2e4, percent_missing=.1, d=13, - ): + +def bgplvm_simulation_missing_data( + optimize=True, + verbose=1, + plot=True, + plot_sim=False, + max_iters=2e4, + percent_missing=0.1, + d=13, +): from GPy import kern from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch @@ -404,28 +491,42 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) - inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype(bool) # 80% missing data + inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype( + bool + ) # 80% missing data Ymissing = Y.copy() Ymissing[inan] = _np.nan - m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, - kernel=k, missing_data=True) + m = BayesianGPLVMMiniBatch( + Ymissing, + Q, + init="random", + num_inducing=num_inducing, + kernel=k, + missing_data=True, + ) m.Yreal = Y if optimize: print("Optimizing model:") - m.optimize('bfgs', messages=verbose, max_iters=max_iters, - gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m -def bgplvm_simulation_missing_data_stochastics(optimize=True, verbose=1, - plot=True, plot_sim=False, - max_iters=2e4, percent_missing=.1, d=13, batchsize=2, - ): + +def bgplvm_simulation_missing_data_stochastics( + optimize=True, + verbose=1, + plot=True, + plot_sim=False, + max_iters=2e4, + percent_missing=0.1, + d=13, + batchsize=2, +): from GPy import kern from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch @@ -434,19 +535,28 @@ def bgplvm_simulation_missing_data_stochastics(optimize=True, verbose=1, Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) - inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype(bool) # 80% missing data + inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype( + bool + ) # 80% missing data Ymissing = Y.copy() Ymissing[inan] = _np.nan - m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, - kernel=k, missing_data=True, stochastic=True, batchsize=batchsize) + m = BayesianGPLVMMiniBatch( + Ymissing, + Q, + init="random", + num_inducing=num_inducing, + kernel=k, + missing_data=True, + stochastic=True, + batchsize=batchsize, + ) m.Yreal = Y if optimize: print("Optimizing model:") - m.optimize('bfgs', messages=verbose, max_iters=max_iters, - gtol=.05) + m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() @@ -461,9 +571,17 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim) k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4) - m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz='permute', **kw) + m = MRD( + Ylist, + input_dim=Q, + num_inducing=num_inducing, + kernel=k, + initx="PCA_concat", + initz="permute", + **kw + ) - m['.*noise'] = [Y.var() / 40. for Y in Ylist] + m[".*noise"] = [Y.var() / 40.0 for Y in Ylist] if optimize: print("Optimizing Model:") @@ -473,7 +591,10 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): m.plot_scales() return m -def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): + +def mrd_simulation_missing_data( + optimize=True, verbose=True, plot=True, plot_sim=True, **kw +): from GPy import kern from GPy.models import MRD @@ -484,29 +605,37 @@ def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim inanlist = [] for Y in Ylist: - inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool) + inan = _np.random.binomial(1, 0.6, size=Y.shape).astype(bool) inanlist.append(inan) Y[inan] = _np.nan - m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, - kernel=k, inference_method=None, - initx="random", initz='permute', **kw) + m = MRD( + Ylist, + input_dim=Q, + num_inducing=num_inducing, + kernel=k, + inference_method=None, + initx="random", + initz="permute", + **kw + ) if optimize: print("Optimizing Model:") - m.optimize('bfgs', messages=verbose, max_iters=8e3, gtol=.1) + m.optimize("bfgs", messages=verbose, max_iters=8e3, gtol=0.1) if plot: m.X.plot("MRD Latent Space 1D") m.plot_scales() return m + def brendan_faces(optimize=True, verbose=True, plot=True): import GPy import pods data = pods.datasets.brendan_faces() Q = 2 - Y = data['Y'] + Y = data["Y"] Yn = Y - Y.mean() Yn /= Yn.std() @@ -514,39 +643,56 @@ def brendan_faces(optimize=True, verbose=True, plot=True): # optimize - if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000) + if optimize: + m.optimize("bfgs", messages=verbose, max_iters=1000) if plot: ax = m.plot_latent(which_indices=(0, 1)) y = m.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False) - lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax) - input('Press enter to finish') + data_show = GPy.plotting.matplot_dep.visualize.image_show( + y[None, :], + dimensions=(20, 28), + transpose=True, + order="F", + invert=False, + scale=False, + ) + lvm = GPy.plotting.matplot_dep.visualize.lvm( + m.X.mean[0, :].copy(), m, data_show, ax + ) + input("Press enter to finish") return m + def olivetti_faces(optimize=True, verbose=True, plot=True): import GPy import pods data = pods.datasets.olivetti_faces() Q = 2 - Y = data['Y'] + Y = data["Y"] Yn = Y - Y.mean() Yn /= Yn.std() m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20) - if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000) + if optimize: + m.optimize("bfgs", messages=verbose, max_iters=1000) if plot: ax = m.plot_latent(which_indices=(0, 1)) y = m.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False) - lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax) - input('Press enter to finish') + data_show = GPy.plotting.matplot_dep.visualize.image_show( + y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False + ) + lvm = GPy.plotting.matplot_dep.visualize.lvm( + m.X.mean[0, :].copy(), m, data_show, ax + ) + input("Press enter to finish") return m + def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=True): import GPy import pods @@ -554,15 +700,18 @@ def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=Tru data = pods.datasets.osu_run1() # optimize if range == None: - Y = data['Y'].copy() + Y = data["Y"].copy() else: - Y = data['Y'][range[0]:range[1], :].copy() + Y = data["Y"][range[0] : range[1], :].copy() if plot: y = Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) + data_show = GPy.plotting.matplot_dep.visualize.stick_show( + y[None, :], connect=data["connect"] + ) GPy.plotting.matplot_dep.visualize.data_play(Y, data_show, frame_rate) return Y + def stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy @@ -570,19 +719,25 @@ def stick(kernel=None, optimize=True, verbose=True, plot=True): data = pods.datasets.osu_run1() # optimize - m = GPy.models.GPLVM(data['Y'], 2, kernel=kernel) - if optimize: m.optimize('bfgs', messages=verbose, max_f_eval=10000) + m = GPy.models.GPLVM(data["Y"], 2, kernel=kernel) + if optimize: + m.optimize("bfgs", messages=verbose, max_f_eval=10000) if plot: plt.clf ax = m.plot_latent() y = m.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) - lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax) - input('Press enter to finish') + data_show = GPy.plotting.matplot_dep.visualize.stick_show( + y[None, :], connect=data["connect"] + ) + lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm( + m.X[:1, :].copy(), m, data_show, latent_axes=ax + ) + input("Press enter to finish") lvm_visualizer.close() data_show.close() return m + def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy @@ -590,19 +745,23 @@ def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True): data = pods.datasets.osu_run1() # optimize - mapping = GPy.mappings.Linear(data['Y'].shape[1], 2) - m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping) - if optimize: m.optimize(messages=verbose, max_f_eval=10000) + mapping = GPy.mappings.Linear(data["Y"].shape[1], 2) + m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping) + if optimize: + m.optimize(messages=verbose, max_f_eval=10000) if plot and GPy.plotting.matplot_dep.visualize.visual_available: plt.clf ax = m.plot_latent() y = m.likelihood.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) + data_show = GPy.plotting.matplot_dep.visualize.stick_show( + y[None, :], connect=data["connect"] + ) GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) - input('Press enter to finish') + input("Press enter to finish") return m + def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy @@ -610,20 +769,24 @@ def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True): data = pods.datasets.osu_run1() # optimize - back_kernel = GPy.kern.RBF(data['Y'].shape[1], lengthscale=5.) - mapping = GPy.mappings.Kernel(X=data['Y'], output_dim=2, kernel=back_kernel) - m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping) - if optimize: m.optimize(messages=verbose, max_f_eval=10000) + back_kernel = GPy.kern.RBF(data["Y"].shape[1], lengthscale=5.0) + mapping = GPy.mappings.Kernel(X=data["Y"], output_dim=2, kernel=back_kernel) + m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping) + if optimize: + m.optimize(messages=verbose, max_f_eval=10000) if plot and GPy.plotting.matplot_dep.visualize.visual_available: plt.clf ax = m.plot_latent() y = m.likelihood.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) + data_show = GPy.plotting.matplot_dep.visualize.stick_show( + y[None, :], connect=data["connect"] + ) GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) # input('Press enter to finish') return m + def robot_wireless(optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy @@ -631,13 +794,15 @@ def robot_wireless(optimize=True, verbose=True, plot=True): data = pods.datasets.robot_wireless() # optimize - m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25) - if optimize: m.optimize(messages=verbose, max_f_eval=10000) + m = GPy.models.BayesianGPLVM(data["Y"], 4, num_inducing=25) + if optimize: + m.optimize(messages=verbose, max_f_eval=10000) if plot: m.plot_latent() return m + def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True): """Interactive visualisation of the Stick Man data from Ohio State University with the Bayesian GPLVM.""" from GPy.models import BayesianGPLVM @@ -648,15 +813,16 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True): data = pods.datasets.osu_run1() Q = 6 - kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True) - m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel) + kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(0.5, Q), ARD=True) + m = BayesianGPLVM(data["Y"], Q, init="PCA", num_inducing=20, kernel=kernel) m.data = data m.likelihood.variance = 0.001 # optimize try: - if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10) + if optimize: + m.optimize("bfgs", messages=verbose, max_iters=5e3, bfgs_factor=10) except KeyboardInterrupt: print("Keyboard interrupt, continuing to plot and return") @@ -665,17 +831,27 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True): plt.sca(latent_axes) m.plot_latent(ax=latent_axes) y = m.Y[:1, :].copy() - data_show = GPy.plotting.matplot_dep.visualize.stick_show(y, connect=data['connect']) - dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes) + data_show = GPy.plotting.matplot_dep.visualize.stick_show( + y, connect=data["connect"] + ) + dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect( + m.X.mean[:1, :].copy(), + m, + data_show, + latent_axes=latent_axes, + sense_axes=sense_axes, + ) fig.canvas.draw() # Canvas.show doesn't work on OSX. - #fig.canvas.show() - input('Press enter to finish') + # fig.canvas.show() + input("Press enter to finish") return m -def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose=True, plot=True): +def cmu_mocap( + subject="35", motion=["01"], in_place=True, optimize=True, verbose=True, plot=True +): import matplotlib.pyplot as plt import GPy import pods @@ -683,34 +859,40 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose data = pods.datasets.cmu_mocap(subject, motion) if in_place: # Make figure move in place. - data['Y'][:, 0:3] = 0.0 - Y = data['Y'] + data["Y"][:, 0:3] = 0.0 + Y = data["Y"] m = GPy.models.GPLVM(Y, 2, normalizer=True) - if optimize: m.optimize(messages=verbose, max_f_eval=10000) + if optimize: + m.optimize(messages=verbose, max_f_eval=10000) if plot: fig, _ = plt.subplots(figsize=(8, 5)) latent_axes = fig.add_subplot(131) sense_axes = fig.add_subplot(132) - viz_axes = fig.add_subplot(133, projection='3d') - + viz_axes = fig.add_subplot(133, projection="3d") m.plot_latent(ax=latent_axes) - latent_axes.set_aspect('equal') + latent_axes.set_aspect("equal") y = m.Y[0, :] - data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel'], viz_axes) + data_show = GPy.plotting.matplot_dep.visualize.skeleton_show( + y[None, :], data["skel"], viz_axes + ) - lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes) - input('Press enter to finish') + lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm( + m.X[0].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes + ) + input("Press enter to finish") lvm_visualizer.close() data_show.close() return m + def ssgplvm_simulation_linear(): import numpy as np import GPy + N, D, Q = 1000, 20, 5 pi = 0.2 @@ -721,7 +903,7 @@ def ssgplvm_simulation_linear(): if dies[q] < pi: x[q] = np.random.randn() else: - x[q] = 0. + x[q] = 0.0 return x Y = np.empty((N, D)) @@ -731,4 +913,3 @@ def ssgplvm_simulation_linear(): X[n] = sample_X(Q, pi) w = np.random.randn(D, Q) Y[n] = np.dot(w, X[n]) - diff --git a/GPy/examples/non_gaussian.py b/GPy/examples/non_gaussian.py index 5d5fc7f6..1f030f2c 100644 --- a/GPy/examples/non_gaussian.py +++ b/GPy/examples/non_gaussian.py @@ -4,38 +4,40 @@ import GPy import numpy as np from GPy.util import datasets + try: import matplotlib.pyplot as plt except: pass + def student_t_approx(optimize=True, plot=True): """ Example of regressing with a student t likelihood using Laplace """ real_std = 0.1 - #Start a function, any function - X = np.linspace(0.0, np.pi*2, 100)[:, None] - Y = np.sin(X) + np.random.randn(*X.shape)*real_std - Y = Y/Y.max() + # Start a function, any function + X = np.linspace(0.0, np.pi * 2, 100)[:, None] + Y = np.sin(X) + np.random.randn(*X.shape) * real_std + Y = Y / Y.max() Yc = Y.copy() - X_full = np.linspace(0.0, np.pi*2, 500)[:, None] + X_full = np.linspace(0.0, np.pi * 2, 500)[:, None] Y_full = np.sin(X_full) - Y_full = Y_full/Y_full.max() + Y_full = Y_full / Y_full.max() - #Slightly noisy data + # Slightly noisy data Yc[75:80] += 1 - #Very noisy data - #Yc[10] += 100 - #Yc[25] += 10 - #Yc[23] += 10 - #Yc[26] += 1000 - #Yc[24] += 10 - #Yc = Yc/Yc.max() + # Very noisy data + # Yc[10] += 100 + # Yc[25] += 10 + # Yc[23] += 10 + # Yc[26] += 1000 + # Yc[24] += 10 + # Yc = Yc/Yc.max() - #Add student t random noise to datapoints + # Add student t random noise to datapoints deg_free = 1 print("Real noise: ", real_std) initial_var_guess = 0.5 @@ -47,45 +49,50 @@ def student_t_approx(optimize=True, plot=True): kernel3 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1]) kernel4 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1]) - #Gaussian GP model on clean data + # Gaussian GP model on clean data m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) # optimize - m1['.*white'].constrain_fixed(1e-5) + m1[".*white"].constrain_fixed(1e-5) m1.randomize() - #Gaussian GP model on corrupt data + # Gaussian GP model on corrupt data m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2) - m2['.*white'].constrain_fixed(1e-5) + m2[".*white"].constrain_fixed(1e-5) m2.randomize() - #Student t GP model on clean data + # Student t GP model on clean data t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd) laplace_inf = GPy.inference.latent_function_inference.Laplace() - m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf) - m3['.*t_scale2'].constrain_bounded(1e-6, 10.) - m3['.*white'].constrain_fixed(1e-5) + m3 = GPy.core.GP( + X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf + ) + m3[".*t_scale2"].constrain_bounded(1e-6, 10.0) + m3[".*white"].constrain_fixed(1e-5) m3.randomize() - #Student t GP model on corrupt data + # Student t GP model on corrupt data t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd) laplace_inf = GPy.inference.latent_function_inference.Laplace() - m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf) - m4['.*t_scale2'].constrain_bounded(1e-6, 10.) - m4['.*white'].constrain_fixed(1e-5) + m4 = GPy.core.GP( + X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf + ) + m4[".*t_scale2"].constrain_bounded(1e-6, 10.0) + m4[".*white"].constrain_fixed(1e-5) m4.randomize() print(m4) - debug=True + debug = True if debug: m4.optimize(messages=1) from matplotlib import pyplot as pb + pb.plot(m4.X, m4.inference_method.f_hat) - pb.plot(m4.X, m4.Y, 'rx') + pb.plot(m4.X, m4.Y, "rx") m4.plot() print(m4) return m4 if optimize: - optimizer='scg' + optimizer = "scg" print("Clean Gaussian") m1.optimize(optimizer, messages=1) print("Corrupt Gaussian") @@ -97,77 +104,91 @@ def student_t_approx(optimize=True, plot=True): if plot: plt.figure(1) - plt.suptitle('Gaussian likelihood') + plt.suptitle("Gaussian likelihood") ax = plt.subplot(211) m1.plot(ax=ax) plt.plot(X_full, Y_full) plt.ylim(-1.5, 1.5) - plt.title('Gaussian clean') + plt.title("Gaussian clean") ax = plt.subplot(212) m2.plot(ax=ax) plt.plot(X_full, Y_full) plt.ylim(-1.5, 1.5) - plt.title('Gaussian corrupt') + plt.title("Gaussian corrupt") plt.figure(2) - plt.suptitle('Student-t likelihood') + plt.suptitle("Student-t likelihood") ax = plt.subplot(211) m3.plot(ax=ax) plt.plot(X_full, Y_full) plt.ylim(-1.5, 1.5) - plt.title('Student-t rasm clean') + plt.title("Student-t rasm clean") ax = plt.subplot(212) m4.plot(ax=ax) plt.plot(X_full, Y_full) plt.ylim(-1.5, 1.5) - plt.title('Student-t rasm corrupt') + plt.title("Student-t rasm corrupt") return m1, m2, m3, m4 + def boston_example(optimize=True, plot=True): raise NotImplementedError("Needs updating") import sklearn from sklearn.cross_validation import KFold - optimizer='bfgs' - messages=0 + + optimizer = "bfgs" + messages = 0 data = datasets.boston_housing() degrees_freedoms = [3, 5, 8, 10] - X = data['X'].copy() - Y = data['Y'].copy() - X = X-X.mean(axis=0) - X = X/X.std(axis=0) - Y = Y-Y.mean() - Y = Y/Y.std() + X = data["X"].copy() + Y = data["Y"].copy() + X = X - X.mean(axis=0) + X = X / X.std(axis=0) + Y = Y - Y.mean() + Y = Y / Y.std() num_folds = 10 kf = KFold(len(Y), n_folds=num_folds, indices=True) - num_models = len(degrees_freedoms) + 3 #3 for baseline, gaussian, gaussian laplace approx + num_models = ( + len(degrees_freedoms) + 3 + ) # 3 for baseline, gaussian, gaussian laplace approx score_folds = np.zeros((num_models, num_folds)) pred_density = score_folds.copy() def rmse(Y, Ystar): - return np.sqrt(np.mean((Y-Ystar)**2)) + return np.sqrt(np.mean((Y - Ystar) ** 2)) for n, (train, test) in enumerate(kf): X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test] print("Fold {}".format(n)) - noise = 1e-1 #np.exp(-2) + noise = 1e-1 # np.exp(-2) rbf_len = 0.5 data_axis_plot = 4 - kernelstu = GPy.kern.RBF(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1]) - kernelgp = GPy.kern.RBF(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1]) + kernelstu = ( + GPy.kern.RBF(X.shape[1]) + + GPy.kern.white(X.shape[1]) + + GPy.kern.bias(X.shape[1]) + ) + kernelgp = ( + GPy.kern.RBF(X.shape[1]) + + GPy.kern.white(X.shape[1]) + + GPy.kern.bias(X.shape[1]) + ) - #Baseline + # Baseline score_folds[0, n] = rmse(Y_test, np.mean(Y_train)) - #Gaussian GP + # Gaussian GP print("Gauss GP") - mgp = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelgp.copy()) - mgp.constrain_fixed('.*white', 1e-5) - mgp['.*len'] = rbf_len - mgp['.*noise'] = noise + mgp = GPy.models.GPRegression( + X_train.copy(), Y_train.copy(), kernel=kernelgp.copy() + ) + mgp.constrain_fixed(".*white", 1e-5) + mgp[".*len"] = rbf_len + mgp[".*noise"] = noise print(mgp) if optimize: mgp.optimize(optimizer=optimizer, messages=messages) @@ -179,13 +200,20 @@ def boston_example(optimize=True, plot=True): print("Gaussian Laplace GP") N, D = Y_train.shape - g_distribution = GPy.likelihoods.noise_model_constructors.gaussian(variance=noise, N=N, D=D) + g_distribution = GPy.likelihoods.noise_model_constructors.gaussian( + variance=noise, N=N, D=D + ) g_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), g_distribution) - mg = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=g_likelihood) - mg.constrain_positive('noise_variance') - mg.constrain_fixed('.*white', 1e-5) - mg['rbf_len'] = rbf_len - mg['noise'] = noise + mg = GPy.models.GPRegression( + X_train.copy(), + Y_train.copy(), + kernel=kernelstu.copy(), + likelihood=g_likelihood, + ) + mg.constrain_positive("noise_variance") + mg.constrain_fixed(".*white", 1e-5) + mg["rbf_len"] = rbf_len + mg["noise"] = noise print(mg) if optimize: mg.optimize(optimizer=optimizer, messages=messages) @@ -196,101 +224,108 @@ def boston_example(optimize=True, plot=True): print(mg) for stu_num, df in enumerate(degrees_freedoms): - #Student T + # Student T print("Student-T GP {}df".format(df)) - t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=df, sigma2=noise) + t_distribution = GPy.likelihoods.noise_model_constructors.student_t( + deg_free=df, sigma2=noise + ) stu_t_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), t_distribution) - mstu_t = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=stu_t_likelihood) - mstu_t.constrain_fixed('.*white', 1e-5) - mstu_t.constrain_bounded('.*t_scale2', 0.0001, 1000) - mstu_t['rbf_len'] = rbf_len - mstu_t['.*t_scale2'] = noise + mstu_t = GPy.models.GPRegression( + X_train.copy(), + Y_train.copy(), + kernel=kernelstu.copy(), + likelihood=stu_t_likelihood, + ) + mstu_t.constrain_fixed(".*white", 1e-5) + mstu_t.constrain_bounded(".*t_scale2", 0.0001, 1000) + mstu_t["rbf_len"] = rbf_len + mstu_t[".*t_scale2"] = noise print(mstu_t) if optimize: mstu_t.optimize(optimizer=optimizer, messages=messages) Y_test_pred = mstu_t.predict(X_test) - score_folds[3+stu_num, n] = rmse(Y_test, Y_test_pred[0]) - pred_density[3+stu_num, n] = np.mean(mstu_t.log_predictive_density(X_test, Y_test)) + score_folds[3 + stu_num, n] = rmse(Y_test, Y_test_pred[0]) + pred_density[3 + stu_num, n] = np.mean( + mstu_t.log_predictive_density(X_test, Y_test) + ) print(pred_density) print(mstu_t) if plot: plt.figure() plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0]) - plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x') - plt.title('GP gauss') + plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x") + plt.title("GP gauss") plt.figure() plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0]) - plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x') - plt.title('Lap gauss') + plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x") + plt.title("Lap gauss") plt.figure() plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0]) - plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x') - plt.title('Stu t {}df'.format(df)) + plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x") + plt.title("Stu t {}df".format(df)) print("Average scores: {}".format(np.mean(score_folds, 1))) print("Average pred density: {}".format(np.mean(pred_density, 1))) if plot: - #Plotting - stu_t_legends = ['Student T, df={}'.format(df) for df in degrees_freedoms] - legends = ['Baseline', 'Gaussian', 'Laplace Approx Gaussian'] + stu_t_legends + # Plotting + stu_t_legends = ["Student T, df={}".format(df) for df in degrees_freedoms] + legends = ["Baseline", "Gaussian", "Laplace Approx Gaussian"] + stu_t_legends - #Plot boxplots for RMSE density + # Plot boxplots for RMSE density fig = plt.figure() - ax=fig.add_subplot(111) - plt.title('RMSE') - bp = ax.boxplot(score_folds.T, notch=0, sym='+', vert=1, whis=1.5) - plt.setp(bp['boxes'], color='black') - plt.setp(bp['whiskers'], color='black') - plt.setp(bp['fliers'], color='red', marker='+') + ax = fig.add_subplot(111) + plt.title("RMSE") + bp = ax.boxplot(score_folds.T, notch=0, sym="+", vert=1, whis=1.5) + plt.setp(bp["boxes"], color="black") + plt.setp(bp["whiskers"], color="black") + plt.setp(bp["fliers"], color="red", marker="+") xtickNames = plt.setp(ax, xticklabels=legends) plt.setp(xtickNames, rotation=45, fontsize=8) - ax.set_ylabel('RMSE') - ax.set_xlabel('Distribution') - #Make grid and put it below boxes - ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', - alpha=0.5) + ax.set_ylabel("RMSE") + ax.set_xlabel("Distribution") + # Make grid and put it below boxes + ax.yaxis.grid(True, linestyle="-", which="major", color="lightgrey", alpha=0.5) ax.set_axisbelow(True) - #Plot boxplots for predictive density + # Plot boxplots for predictive density fig = plt.figure() - ax=fig.add_subplot(111) - plt.title('Predictive density') - bp = ax.boxplot(pred_density[1:,:].T, notch=0, sym='+', vert=1, whis=1.5) - plt.setp(bp['boxes'], color='black') - plt.setp(bp['whiskers'], color='black') - plt.setp(bp['fliers'], color='red', marker='+') + ax = fig.add_subplot(111) + plt.title("Predictive density") + bp = ax.boxplot(pred_density[1:, :].T, notch=0, sym="+", vert=1, whis=1.5) + plt.setp(bp["boxes"], color="black") + plt.setp(bp["whiskers"], color="black") + plt.setp(bp["fliers"], color="red", marker="+") xtickNames = plt.setp(ax, xticklabels=legends[1:]) plt.setp(xtickNames, rotation=45, fontsize=8) - ax.set_ylabel('Mean Log probability P(Y*|Y)') - ax.set_xlabel('Distribution') - #Make grid and put it below boxes - ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey', - alpha=0.5) + ax.set_ylabel("Mean Log probability P(Y*|Y)") + ax.set_xlabel("Distribution") + # Make grid and put it below boxes + ax.yaxis.grid(True, linestyle="-", which="major", color="lightgrey", alpha=0.5) ax.set_axisbelow(True) return mstu_t -#def precipitation_example(): - #import sklearn - #from sklearn.cross_validation import KFold - #data = datasets.boston_housing() - #X = data['X'].copy() - #Y = data['Y'].copy() - #X = X-X.mean(axis=0) - #X = X/X.std(axis=0) - #Y = Y-Y.mean() - #Y = Y/Y.std() - #import ipdb; ipdb.set_trace() # XXX BREAKPOINT - #num_folds = 10 - #kf = KFold(len(Y), n_folds=num_folds, indices=True) - #score_folds = np.zeros((4, num_folds)) - #def rmse(Y, Ystar): - #return np.sqrt(np.mean((Y-Ystar)**2)) - ##for train, test in kf: - #for n, (train, test) in enumerate(kf): - #X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test] - #print "Fold {}".format(n) +# def precipitation_example(): +# import sklearn +# from sklearn.cross_validation import KFold +# data = datasets.boston_housing() +# X = data['X'].copy() +# Y = data['Y'].copy() +# X = X-X.mean(axis=0) +# X = X/X.std(axis=0) +# Y = Y-Y.mean() +# Y = Y/Y.std() +# import ipdb; ipdb.set_trace() # XXX BREAKPOINT +# num_folds = 10 +# kf = KFold(len(Y), n_folds=num_folds, indices=True) +# score_folds = np.zeros((4, num_folds)) +# def rmse(Y, Ystar): +# return np.sqrt(np.mean((Y-Ystar)**2)) +##for train, test in kf: +# for n, (train, test) in enumerate(kf): +# X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test] +# print "Fold {}".format(n) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 93e05cbb..9ade6552 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -11,88 +11,112 @@ except: import numpy as np import GPy + def olympic_marathon_men(optimize=True, plot=True): """Run a standard Gaussian process regression on the Olympic marathon data.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.olympic_marathon_men() # create simple GP Model - m = GPy.models.GPRegression(data['X'], data['Y']) + m = GPy.models.GPRegression(data["X"], data["Y"]) # set the lengthscale to be something sensible (defaults to 1) - m.kern.lengthscale = 10. + m.kern.lengthscale = 10.0 if optimize: - m.optimize('bfgs', max_iters=200) + m.optimize("bfgs", max_iters=200) if plot: m.plot(plot_limits=(1850, 2050)) return m + def coregionalization_toy(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions. """ - #build a design matrix with a column of integers indicating the output + # build a design matrix with a column of integers indicating the output X1 = np.random.rand(50, 1) * 8 X2 = np.random.rand(30, 1) * 5 - #build a suitable set of observed variables + # build a suitable set of observed variables Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 - Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. + Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.0 - m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) + m = GPy.models.GPCoregionalizedRegression(X_list=[X1, X2], Y_list=[Y1, Y2]) if optimize: - m.optimize('bfgs', max_iters=100) + m.optimize("bfgs", max_iters=100) if plot: - slices = GPy.util.multioutput.get_slices([X1,X2]) - m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) - m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) + slices = GPy.util.multioutput.get_slices([X1, X2]) + m.plot( + fixed_inputs=[(1, 0)], + which_data_rows=slices[0], + Y_metadata={"output_index": 0}, + ) + m.plot( + fixed_inputs=[(1, 1)], + which_data_rows=slices[1], + Y_metadata={"output_index": 1}, + ax=pb.gca(), + ) return m + def coregionalization_sparse(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations. """ - #build a design matrix with a column of integers indicating the output + # build a design matrix with a column of integers indicating the output X1 = np.random.rand(50, 1) * 8 X2 = np.random.rand(30, 1) * 5 - #build a suitable set of observed variables + # build a suitable set of observed variables Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 - Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. + Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.0 - m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) + m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1, X2], Y_list=[Y1, Y2]) if optimize: - m.optimize('bfgs', max_iters=100) + m.optimize("bfgs", max_iters=100) if plot: - slices = GPy.util.multioutput.get_slices([X1,X2]) - m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) - m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) + slices = GPy.util.multioutput.get_slices([X1, X2]) + m.plot( + fixed_inputs=[(1, 0)], + which_data_rows=slices[0], + Y_metadata={"output_index": 0}, + ) + m.plot( + fixed_inputs=[(1, 1)], + which_data_rows=slices[1], + Y_metadata={"output_index": 1}, + ax=pb.gca(), + ) pb.ylim(-3,) return m + def epomeo_gpx(max_iters=200, optimize=True, plot=True): """ Perform Gaussian process regression on the latitude and longitude data from the Mount Epomeo runs. Requires gpxpy to be installed on your system to load in the data. """ - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.epomeo_gpx() num_data_list = [] - for Xpart in data['X']: + for Xpart in data["X"]: num_data_list.append(Xpart.shape[0]) num_data_array = np.array(num_data_list) @@ -100,29 +124,43 @@ def epomeo_gpx(max_iters=200, optimize=True, plot=True): Y = np.zeros((num_data, 2)) t = np.zeros((num_data, 2)) start = 0 - for Xpart, index in zip(data['X'], range(len(data['X']))): - end = start+Xpart.shape[0] - t[start:end, :] = np.hstack((Xpart[:, 0:1], - index*np.ones((Xpart.shape[0], 1)))) + for Xpart, index in zip(data["X"], range(len(data["X"]))): + end = start + Xpart.shape[0] + t[start:end, :] = np.hstack( + (Xpart[:, 0:1], index * np.ones((Xpart.shape[0], 1))) + ) Y[start:end, :] = Xpart[:, 1:3] num_inducing = 200 - Z = np.hstack((np.linspace(t[:,0].min(), t[:, 0].max(), num_inducing)[:, None], - np.random.randint(0, 4, num_inducing)[:, None])) + Z = np.hstack( + ( + np.linspace(t[:, 0].min(), t[:, 0].max(), num_inducing)[:, None], + np.random.randint(0, 4, num_inducing)[:, None], + ) + ) k1 = GPy.kern.RBF(1) k2 = GPy.kern.Coregionalize(output_dim=5, rank=5) - k = k1**k2 + k = k1 ** k2 m = GPy.models.SparseGPRegression(t, Y, kernel=k, Z=Z, normalize_Y=True) - m.constrain_fixed('.*variance', 1.) + m.constrain_fixed(".*variance", 1.0) m.inducing_inputs.constrain_fixed() m.Gaussian_noise.variance.constrain_bounded(1e-3, 1e-1) - m.optimize(max_iters=max_iters,messages=True) + m.optimize(max_iters=max_iters, messages=True) return m -def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True): + +def multiple_optima( + gene_number=937, + resolution=80, + model_restarts=10, + seed=10000, + max_iters=300, + optimize=True, + plot=True, +): """ Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisy mode is @@ -130,25 +168,30 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000 """ # Contour over a range of length scales and signal/noise ratios. - length_scales = np.linspace(0.1, 60., resolution) - log_SNRs = np.linspace(-3., 4., resolution) + length_scales = np.linspace(0.1, 60.0, resolution) + log_SNRs = np.linspace(-3.0, 4.0, resolution) - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return - data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number) + data = pods.datasets.della_gatta_TRP63_gene_expression( + data_set="della_gatta", gene_number=gene_number + ) # data['Y'] = data['Y'][0::2, :] # data['X'] = data['X'][0::2, :] - data['Y'] = data['Y'] - np.mean(data['Y']) + data["Y"] = data["Y"] - np.mean(data["Y"]) - lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF) + lls = GPy.examples.regression._contour_data( + data, length_scales, log_SNRs, GPy.kern.RBF + ) if plot: pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet) ax = pb.gca() - pb.xlabel('length scale') - pb.ylabel('log_10 SNR') + pb.xlabel("length scale") + pb.ylabel("log_10 SNR") xlim = ax.get_xlim() ylim = ax.get_ylim() @@ -160,28 +203,41 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000 np.random.seed(seed=seed) for i in range(0, model_restarts): # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) - kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50)) + kern = GPy.kern.RBF( + 1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50) + ) - m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern) + m = GPy.models.GPRegression(data["X"], data["Y"], kernel=kern) m.likelihood.variance = np.random.uniform(1e-3, 1) optim_point_x[0] = m.rbf.lengthscale - optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance); + optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance) # optimize if optimize: - m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters) + m.optimize("scg", xtol=1e-6, ftol=1e-6, max_iters=max_iters) optim_point_x[1] = m.rbf.lengthscale - optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance); + optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance) if plot: - pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k') + pb.arrow( + optim_point_x[0], + optim_point_y[0], + optim_point_x[1] - optim_point_x[0], + optim_point_y[1] - optim_point_y[0], + label=str(i), + head_length=1, + head_width=0.5, + fc="k", + ec="k", + ) models.append(m) if plot: ax.set_xlim(xlim) ax.set_ylim(ylim) - return m # (models, lls) + return m # (models, lls) + def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF): """ @@ -195,19 +251,19 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF): """ lls = [] - total_var = np.var(data['Y']) - kernel = kernel_call(1, variance=1., lengthscale=1.) - model = GPy.models.GPRegression(data['X'], data['Y'], kernel=kernel) + total_var = np.var(data["Y"]) + kernel = kernel_call(1, variance=1.0, lengthscale=1.0) + model = GPy.models.GPRegression(data["X"], data["Y"], kernel=kernel) for log_SNR in log_SNRs: - SNR = 10.**log_SNR - noise_var = total_var / (1. + SNR) + SNR = 10.0 ** log_SNR + noise_var = total_var / (1.0 + SNR) signal_var = total_var - noise_var - model.kern['.*variance'] = signal_var + model.kern[".*variance"] = signal_var model.likelihood.variance = noise_var length_scale_lls = [] for length_scale in length_scales: - model['.*lengthscale'] = length_scale + model[".*lengthscale"] = length_scale length_scale_lls.append(model.log_likelihood()) lls.append(length_scale_lls) @@ -217,86 +273,97 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF): def olympic_100m_men(optimize=True, plot=True): """Run a standard Gaussian process regression on the Rogers and Girolami olympics data.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.olympic_100m_men() # create simple GP Model - m = GPy.models.GPRegression(data['X'], data['Y']) + m = GPy.models.GPRegression(data["X"], data["Y"]) # set the lengthscale to be something sensible (defaults to 1) m.rbf.lengthscale = 10 if optimize: - m.optimize('bfgs', max_iters=200) + m.optimize("bfgs", max_iters=200) if plot: m.plot(plot_limits=(1850, 2050)) return m + def toy_rbf_1d(optimize=True, plot=True): """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.toy_rbf_1d() # create simple GP Model - m = GPy.models.GPRegression(data['X'], data['Y']) + m = GPy.models.GPRegression(data["X"], data["Y"]) if optimize: - m.optimize('bfgs') + m.optimize("bfgs") if plot: m.plot() return m + def toy_rbf_1d_50(optimize=True, plot=True): """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.toy_rbf_1d_50() # create simple GP Model - m = GPy.models.GPRegression(data['X'], data['Y']) + m = GPy.models.GPRegression(data["X"], data["Y"]) if optimize: - m.optimize('bfgs') + m.optimize("bfgs") if plot: m.plot() return m + def toy_poisson_rbf_1d_laplace(optimize=True, plot=True): """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" - optimizer='scg' + optimizer = "scg" x_len = 100 X = np.linspace(0, 10, x_len)[:, None] f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.RBF(1).K(X)) - Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None] + Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:, None] kern = GPy.kern.RBF(1) poisson_lik = GPy.likelihoods.Poisson() laplace_inf = GPy.inference.latent_function_inference.Laplace() # create simple GP Model - m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf) + m = GPy.core.GP( + X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf + ) if optimize: m.optimize(optimizer) if plot: m.plot() # plot the real underlying rate function - pb.plot(X, np.exp(f_true), '--k', linewidth=2) + pb.plot(X, np.exp(f_true), "--k", linewidth=2) return m -def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True): + +def toy_ARD( + max_iters=1000, kernel_type="linear", num_samples=300, D=4, optimize=True, plot=True +): # Create an artificial dataset where the values in the targets (Y) # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to # see if this dependency can be recovered @@ -310,14 +377,14 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0])).reshape(-1, 1) Y = np.hstack((Y1, Y2)) - Y = np.dot(Y, np.random.rand(2, D)); + Y = np.dot(Y, np.random.rand(2, D)) Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1]) Y -= Y.mean() Y /= Y.std() - if kernel_type == 'linear': + if kernel_type == "linear": kernel = GPy.kern.Linear(X.shape[1], ARD=1) - elif kernel_type == 'rbf_inv': + elif kernel_type == "rbf_inv": kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1) else: kernel = GPy.kern.RBF(X.shape[1], ARD=1) @@ -327,14 +394,17 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize # m.set_prior('.*lengthscale',len_prior) if optimize: - m.optimize(optimizer='scg', max_iters=max_iters) + m.optimize(optimizer="scg", max_iters=max_iters) if plot: m.kern.plot_ARD() return m -def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True): + +def toy_ARD_sparse( + max_iters=1000, kernel_type="linear", num_samples=300, D=4, optimize=True, plot=True +): # Create an artificial dataset where the values in the targets (Y) # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to # see if this dependency can be recovered @@ -348,69 +418,73 @@ def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, o Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0]))[:, None] Y = np.hstack((Y1, Y2)) - Y = np.dot(Y, np.random.rand(2, D)); + Y = np.dot(Y, np.random.rand(2, D)) Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1]) Y -= Y.mean() Y /= Y.std() - if kernel_type == 'linear': + if kernel_type == "linear": kernel = GPy.kern.Linear(X.shape[1], ARD=1) - elif kernel_type == 'rbf_inv': + elif kernel_type == "rbf_inv": kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1) else: kernel = GPy.kern.RBF(X.shape[1], ARD=1) - #kernel += GPy.kern.Bias(X.shape[1]) + # kernel += GPy.kern.Bias(X.shape[1]) X_variance = np.ones(X.shape) * 0.5 m = GPy.models.SparseGPRegression(X, Y, kernel, X_variance=X_variance) # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25 # m.set_prior('.*lengthscale',len_prior) if optimize: - m.optimize(optimizer='scg', max_iters=max_iters) + m.optimize(optimizer="scg", max_iters=max_iters) if plot: m.kern.plot_ARD() return m + def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True): """Predict the location of a robot given wirelss signal strength readings.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.robot_wireless() # create simple GP Model - m = GPy.models.GPRegression(data['Y'], data['X'], kernel=kernel) + m = GPy.models.GPRegression(data["Y"], data["X"], kernel=kernel) # optimize if optimize: m.optimize(max_iters=max_iters) - Xpredict = m.predict(data['Ytest'])[0] + Xpredict = m.predict(data["Ytest"])[0] if plot: - pb.plot(data['Xtest'][:, 0], data['Xtest'][:, 1], 'r-') - pb.plot(Xpredict[:, 0], Xpredict[:, 1], 'b-') - pb.axis('equal') - pb.title('WiFi Localization with Gaussian Processes') - pb.legend(('True Location', 'Predicted Location')) + pb.plot(data["Xtest"][:, 0], data["Xtest"][:, 1], "r-") + pb.plot(Xpredict[:, 0], Xpredict[:, 1], "b-") + pb.axis("equal") + pb.title("WiFi Localization with Gaussian Processes") + pb.legend(("True Location", "Predicted Location")) - sse = ((data['Xtest'] - Xpredict)**2).sum() + sse = ((data["Xtest"] - Xpredict) ** 2).sum() - print(('Sum of squares error on test data: ' + str(sse))) + print(("Sum of squares error on test data: " + str(sse))) return m + def silhouette(max_iters=100, optimize=True, plot=True): """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper.""" - try:import pods + try: + import pods except ImportError: - print('pods unavailable, see https://github.com/sods/ods for example datasets') + print("pods unavailable, see https://github.com/sods/ods for example datasets") return data = pods.datasets.silhouette() # create simple GP Model - m = GPy.models.GPRegression(data['X'], data['Y']) + m = GPy.models.GPRegression(data["X"], data["Y"]) # optimize if optimize: @@ -419,10 +493,18 @@ def silhouette(max_iters=100, optimize=True, plot=True): print(m) return m -def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, optimize=True, plot=True, checkgrad=False): + +def sparse_GP_regression_1D( + num_samples=400, + num_inducing=5, + max_iters=100, + optimize=True, + plot=True, + checkgrad=False, +): """Run a 1D example of a sparse GP regression.""" # sample inputs and outputs - X = np.random.uniform(-3., 3., (num_samples, 1)) + X = np.random.uniform(-3.0, 3.0, (num_samples, 1)) Y = np.sin(X) + np.random.randn(num_samples, 1) * 0.05 # construct kernel rbf = GPy.kern.RBF(1) @@ -433,20 +515,23 @@ def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, opti m.checkgrad() if optimize: - m.optimize('tnc', max_iters=max_iters) + m.optimize("tnc", max_iters=max_iters) if plot: m.plot() return m -def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False): + +def sparse_GP_regression_2D( + num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False +): """Run a 2D example of a sparse GP regression.""" np.random.seed(1234) - X = np.random.uniform(-3., 3., (num_samples, 2)) + X = np.random.uniform(-3.0, 3.0, (num_samples, 2)) Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05 if nan: - inan = np.random.binomial(1,.2,size=Y.shape) + inan = np.random.binomial(1, 0.2, size=Y.shape) Y[inan] = np.nan # construct kernel @@ -456,13 +541,13 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing) # contrain all parameters to be positive (but not inducing inputs) - m['.*len'] = 2. + m[".*len"] = 2.0 m.checkgrad() # optimize if optimize: - m.optimize('tnc', messages=1, max_iters=max_iters) + m.optimize("tnc", messages=1, max_iters=max_iters) # plot if plot: @@ -471,76 +556,79 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt print(m) return m + def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): """Run a 1D example of a sparse GP regression with uncertain inputs.""" fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True) # sample inputs and outputs S = np.ones((20, 1)) - X = np.random.uniform(-3., 3., (20, 1)) + X = np.random.uniform(-3.0, 3.0, (20, 1)) Y = np.sin(X) + np.random.randn(20, 1) * 0.05 # likelihood = GPy.likelihoods.Gaussian(Y) - Z = np.random.uniform(-3., 3., (7, 1)) + Z = np.random.uniform(-3.0, 3.0, (7, 1)) k = GPy.kern.RBF(1) # create simple GP Model - no input uncertainty on this one m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z) if optimize: - m.optimize('scg', messages=1, max_iters=max_iters) + m.optimize("scg", messages=1, max_iters=max_iters) if plot: m.plot(ax=axes[0]) - axes[0].set_title('no input uncertainty') + axes[0].set_title("no input uncertainty") print(m) # the same Model with uncertainty m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z, X_variance=S) if optimize: - m.optimize('scg', messages=1, max_iters=max_iters) + m.optimize("scg", messages=1, max_iters=max_iters) if plot: m.plot(ax=axes[1]) - axes[1].set_title('with input uncertainty') + axes[1].set_title("with input uncertainty") fig.canvas.draw() print(m) return m + def simple_mean_function(max_iters=100, optimize=True, plot=True): """ The simplest possible mean function. No parameters, just a simple Sinusoid. """ - #create simple mean function - mf = GPy.core.Mapping(1,1) + # create simple mean function + mf = GPy.core.Mapping(1, 1) mf.f = np.sin - mf.update_gradients = lambda a,b: None + mf.update_gradients = lambda a, b: None - X = np.linspace(0,10,50).reshape(-1,1) - Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) + X = np.linspace(0, 10, 50).reshape(-1, 1) + Y = np.sin(X) + 0.5 * np.cos(3 * X) + 0.1 * np.random.randn(*X.shape) - k =GPy.kern.RBF(1) + k = GPy.kern.RBF(1) lik = GPy.likelihoods.Gaussian() m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf) if optimize: m.optimize(max_iters=max_iters) if plot: - m.plot(plot_limits=(-10,15)) + m.plot(plot_limits=(-10, 15)) return m + def parametric_mean_function(max_iters=100, optimize=True, plot=True): """ A linear mean function with parameters that we'll learn alongside the kernel """ - #create simple mean function - mf = GPy.core.Mapping(1,1) + # create simple mean function + mf = GPy.core.Mapping(1, 1) mf.f = np.sin - X = np.linspace(0,10,50).reshape(-1,1) - Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) + 3*X + X = np.linspace(0, 10, 50).reshape(-1, 1) + Y = np.sin(X) + 0.5 * np.cos(3 * X) + 0.1 * np.random.randn(*X.shape) + 3 * X - mf = GPy.mappings.Linear(1,1) + mf = GPy.mappings.Linear(1, 1) - k =GPy.kern.RBF(1) + k = GPy.kern.RBF(1) lik = GPy.likelihoods.Gaussian() m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf) if optimize: @@ -556,23 +644,27 @@ def warped_gp_cubic_sine(max_iters=100): Snelson's paper. """ X = (2 * np.pi) * np.random.random(151) - np.pi - Y = np.sin(X) + np.random.normal(0,0.2,151) - Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y]) + Y = np.sin(X) + np.random.normal(0, 0.2, 151) + Y = np.array([np.power(abs(y), float(1) / 3) * (1, -1)[y < 0] for y in Y]) X = X[:, None] Y = Y[:, None] warp_k = GPy.kern.RBF(1) warp_f = GPy.util.warping_functions.TanhFunction(n_terms=2) warp_m = GPy.models.WarpedGP(X, Y, kernel=warp_k, warping_function=warp_f) - warp_m['.*\.d'].constrain_fixed(1.0) + warp_m[".*\.d"].constrain_fixed(1.0) m = GPy.models.GPRegression(X, Y) - m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters) - warp_m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters) - #m.optimize(max_iters=max_iters) - #warp_m.optimize(max_iters=max_iters) + m.optimize_restarts( + parallel=False, robust=True, num_restarts=5, max_iters=max_iters + ) + warp_m.optimize_restarts( + parallel=False, robust=True, num_restarts=5, max_iters=max_iters + ) + # m.optimize(max_iters=max_iters) + # warp_m.optimize(max_iters=max_iters) print(warp_m) - print(warp_m['.*warp.*']) + print(warp_m[".*warp.*"]) warp_m.predict_in_warped_space = False warp_m.plot(title="Warped GP - Latent space") @@ -584,55 +676,88 @@ def warped_gp_cubic_sine(max_iters=100): return warp_m - def multioutput_gp_with_derivative_observations(): - def plot_gp_vs_real(m, x, yreal, size_inputs, title, fixed_input=1, xlim=[0,11], ylim=[-1.5,3]): + def plot_gp_vs_real( + m, x, yreal, size_inputs, title, fixed_input=1, xlim=[0, 11], ylim=[-1.5, 3] + ): fig, ax = pb.subplots() ax.set_title(title) - pb.plot(x, yreal, "r", label='Real function') - rows = slice(0, size_inputs[0]) if fixed_input == 0 else slice(size_inputs[0], size_inputs[0]+size_inputs[1]) - m.plot(fixed_inputs=[(1, fixed_input)], which_data_rows=rows, xlim=xlim, ylim=ylim, ax=ax) - f = lambda x: np.sin(x)+0.1*(x-2.)**2-0.005*x**3 - fd = lambda x: np.cos(x)+0.2*(x-2.)-0.015*x**2 - N=10 # Number of observations - M=10 # Number of derivative observations - Npred=100 # Number of prediction points - sigma = 0.05 # Noise of observations - sigma_der = 0.05 # Noise of derivative observations - x = np.array([np.linspace(1,10,N)]).T - y = f(x) + np.array(sigma*np.random.normal(0,1,(N,1))) + pb.plot(x, yreal, "r", label="Real function") + rows = ( + slice(0, size_inputs[0]) + if fixed_input == 0 + else slice(size_inputs[0], size_inputs[0] + size_inputs[1]) + ) + m.plot( + fixed_inputs=[(1, fixed_input)], + which_data_rows=rows, + xlim=xlim, + ylim=ylim, + ax=ax, + ) - xd = np.array([np.linspace(2,8,M)]).T - yd = fd(xd) + np.array(sigma_der*np.random.normal(0,1,(M,1))) + f = lambda x: np.sin(x) + 0.1 * (x - 2.0) ** 2 - 0.005 * x ** 3 + fd = lambda x: np.cos(x) + 0.2 * (x - 2.0) - 0.015 * x ** 2 + N = 10 # Number of observations + M = 10 # Number of derivative observations + Npred = 100 # Number of prediction points + sigma = 0.05 # Noise of observations + sigma_der = 0.05 # Noise of derivative observations + x = np.array([np.linspace(1, 10, N)]).T + y = f(x) + np.array(sigma * np.random.normal(0, 1, (N, 1))) - xpred = np.array([np.linspace(0,11,Npred)]).T + xd = np.array([np.linspace(2, 8, M)]).T + yd = fd(xd) + np.array(sigma_der * np.random.normal(0, 1, (M, 1))) + + xpred = np.array([np.linspace(0, 11, Npred)]).T ypred_true = f(xpred) ydpred_true = fd(xpred) # squared exponential kernel: - se = GPy.kern.RBF(input_dim = 1, lengthscale=1.5, variance=0.2) + se = GPy.kern.RBF(input_dim=1, lengthscale=1.5, variance=0.2) # We need to generate separate kernel for the derivative observations and give the created kernel as an input: se_der = GPy.kern.DiffKern(se, 0) - #Then - gauss = GPy.likelihoods.Gaussian(variance=sigma**2) - gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der**2) + # Then + gauss = GPy.likelihoods.Gaussian(variance=sigma ** 2) + gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der ** 2) # Then create the model, we give everything in lists, the order of the inputs indicates the order of the outputs # Now we have the regular observations first and derivative observations second, meaning that the kernels and # the likelihoods must follow the same order. Crosscovariances are automatically taken care of - m = GPy.models.MultioutputGP(X_list=[x, xd], Y_list=[y, yd], - kernel_list=[se, se_der], - likelihood_list=[gauss, gauss_der]) + m = GPy.models.MultioutputGP( + X_list=[x, xd], + Y_list=[y, yd], + kernel_list=[se, se_der], + likelihood_list=[gauss, gauss_der], + ) # Optimize the model m.optimize(messages=0, ipython_notebook=False) - #Plot the model, the syntax is same as for multioutput models: - plot_gp_vs_real(m, xpred, ydpred_true, [x.shape[0], xd.shape[0]], title='Latent function derivatives', fixed_input=1, xlim=[0,11], ylim=[-1.5,3]) - plot_gp_vs_real(m, xpred, ypred_true, [x.shape[0], xd.shape[0]], title='Latent function', fixed_input=0, xlim=[0,11], ylim=[-1.5,3]) + # Plot the model, the syntax is same as for multioutput models: + plot_gp_vs_real( + m, + xpred, + ydpred_true, + [x.shape[0], xd.shape[0]], + title="Latent function derivatives", + fixed_input=1, + xlim=[0, 11], + ylim=[-1.5, 3], + ) + plot_gp_vs_real( + m, + xpred, + ypred_true, + [x.shape[0], xd.shape[0]], + title="Latent function", + fixed_input=0, + xlim=[0, 11], + ylim=[-1.5, 3], + ) - #making predictions for the values: - mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0,1))]) + # making predictions for the values: + mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0, 1))]) return m diff --git a/GPy/examples/state_space.py b/GPy/examples/state_space.py index 898d1676..33faf89d 100644 --- a/GPy/examples/state_space.py +++ b/GPy/examples/state_space.py @@ -4,26 +4,26 @@ import matplotlib.pyplot as plt import GPy.models.state_space_model as SS_model + def state_space_example(): X = np.linspace(0, 10, 2000)[:, None] - Y = np.sin(X) + np.random.randn(*X.shape)*0.1 + Y = np.sin(X) + np.random.randn(*X.shape) * 0.1 kernel1 = GPy.kern.Matern32(X.shape[1]) - m1 = GPy.models.GPRegression(X,Y, kernel1) + m1 = GPy.models.GPRegression(X, Y, kernel1) print(m1) - m1.optimize(optimizer='bfgs',messages=True) + m1.optimize(optimizer="bfgs", messages=True) print(m1) kernel2 = GPy.kern.sde_Matern32(X.shape[1]) - #m2 = SS_model.StateSpace(X,Y, kernel2) - m2 = GPy.models.StateSpace(X,Y, kernel2) + # m2 = SS_model.StateSpace(X,Y, kernel2) + m2 = GPy.models.StateSpace(X, Y, kernel2) print(m2) - m2.optimize(optimizer='bfgs',messages=True) + m2.optimize(optimizer="bfgs", messages=True) print(m2) return m1, m2 -