From bef114eabd592fc2b44bb823c0e8ec779f25306d Mon Sep 17 00:00:00 2001 From: beckdaniel Date: Tue, 8 Dec 2015 13:59:46 +0000 Subject: [PATCH] (wpgs) fixing newton-raphson for f_inv and fixing plotting stuff --- GPy/models/warped_gp.py | 37 +++++++++++++++++++++--------- GPy/plotting/gpy_plot/plot_util.py | 8 ++++++- GPy/testing/model_tests.py | 29 ++++++++++------------- GPy/util/warping_functions.py | 6 ++++- 4 files changed, 50 insertions(+), 30 deletions(-) diff --git a/GPy/models/warped_gp.py b/GPy/models/warped_gp.py index df947d3e..998c0ed6 100644 --- a/GPy/models/warped_gp.py +++ b/GPy/models/warped_gp.py @@ -96,11 +96,14 @@ class WarpedGP(GP): return arg1 - (arg2 ** 2) def predict(self, Xnew, which_parts='all', pred_init=None, full_cov=False, Y_metadata=None, - median=False, deg_gauss_hermite=100): - # normalize X values - # Xnew = (Xnew.copy() - self._Xoffset) / self._Xscale + median=False, deg_gauss_hermite=100, likelihood=None): + """ + Prediction results depend on: + - The value of the self.predict_in_warped_space flag + - The median flag passed as argument + The likelihood keyword is never used, it is just to follow the plotting API. + """ mu, var = GP._raw_predict(self, Xnew) - # now push through likelihood mean, var = self.likelihood.predictive_values(mu, var) @@ -116,13 +119,11 @@ class WarpedGP(GP): else: wmean = mean wvar = var - if self.scale_data: pred = self._unscale_data(pred) - return wmean, wvar - def predict_quantiles(self, X, quantiles=(2.5, 97.5), Y_metadata=None): + def predict_quantiles(self, X, quantiles=(2.5, 97.5), Y_metadata=None, likelihood=None, median=False): """ Get the predictive quantiles around the prediction at X @@ -137,15 +138,29 @@ class WarpedGP(GP): if self.normalizer is not None: m, v = self.normalizer.inverse_mean(m), self.normalizer.inverse_variance(v) a, b = self.likelihood.predictive_quantiles(m, v, quantiles, Y_metadata) - #return [a, b] if not self.predict_in_warped_space: return [a, b] - #print a.shape new_a = self.warping_function.f_inv(a) new_b = self.warping_function.f_inv(b) - return [new_a, new_b] - #return self.likelihood.predictive_quantiles(m, v, quantiles, Y_metadata) + + def log_predictive_density(self, x_test, y_test, Y_metadata=None): + """ + Calculation of the log predictive density + + .. math: + p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*}) + + :param x_test: test locations (x_{*}) + :type x_test: (Nx1) array + :param y_test: test observations (y_{*}) + :type y_test: (Nx1) array + :param Y_metadata: metadata associated with the test points + """ + mu_star, var_star = self._raw_predict(x_test) + fy = self.warping_function.f(y_test) + ll_lpd = self.likelihood.log_predictive_density(fy, mu_star, var_star, Y_metadata=Y_metadata) + return ll_lpd * self.warping_function.fgrad_y(y_test) if __name__ == '__main__': diff --git a/GPy/plotting/gpy_plot/plot_util.py b/GPy/plotting/gpy_plot/plot_util.py index 254886a2..8c5c04f1 100644 --- a/GPy/plotting/gpy_plot/plot_util.py +++ b/GPy/plotting/gpy_plot/plot_util.py @@ -31,6 +31,7 @@ import numpy as np from scipy import sparse import itertools +from GPy.models import WarpedGP def in_ipynb(): try: @@ -73,6 +74,9 @@ def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, wh if 'output_index' not in predict_kw['Y_metadata']: predict_kw['Y_metadata']['output_index'] = Xgrid[:,-1:].astype(np.int) + if isinstance(self, WarpedGP) and self.predict_in_warped_space: + predict_kw['median'] = True + mu, _ = self.predict(Xgrid, **predict_kw) if percentiles is not None: @@ -291,6 +295,8 @@ def get_x_y_var(model): Y = model.Y.values except AttributeError: Y = model.Y + if isinstance(model, WarpedGP) and model.predict_in_warped_space: + Y = model.Y_untransformed if sparse.issparse(Y): Y = Y.todense().view(np.ndarray) return X, X_variance, Y @@ -377,4 +383,4 @@ def x_frame2D(X,plot_limits=None,resolution=None): resolution = resolution or 50 xx, yy = np.mgrid[xmin[0]:xmax[0]:1j*resolution,xmin[1]:xmax[1]:1j*resolution] Xnew = np.vstack((xx.flatten(),yy.flatten())).T - return Xnew, xx, yy, xmin, xmax \ No newline at end of file + return Xnew, xx, yy, xmin, xmax diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py index 1212d746..988bde13 100644 --- a/GPy/testing/model_tests.py +++ b/GPy/testing/model_tests.py @@ -307,42 +307,37 @@ class MiscTests(unittest.TestCase): np.testing.assert_almost_equal(preds, warp_preds) - @unittest.skip('Comment this to plot the modified sine function') + #@unittest.skip('Comment this to plot the modified sine function') def test_warped_gp_sine(self): """ A test replicating the sine regression problem from Snelson's paper. """ X = (2 * np.pi) * np.random.random(151) - np.pi - Y = np.sin(X) + np.random.normal(0,0.1,151) - Y = np.exp(Y) - 5 - #Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y]) + 0 + 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]) - #np.seterr(over='raise') import matplotlib.pyplot as plt warp_k = GPy.kern.RBF(1) warp_f = GPy.util.warping_functions.TanhWarpingFunction_d(n_terms=2) warp_m = GPy.models.WarpedGP(X[:, None], Y[:, None], kernel=warp_k, warping_function=warp_f) - #warp_m['.*variance.*'].constrain_fixed(0.25) - #warp_m['.*lengthscale.*'].constrain_fixed(1) - #warp_m['warp_tanh.d'].constrain_fixed(1) - #warp_m.randomize() - #warp_m['.*warp_tanh.psi*'][:,0:2].constrain_bounded(0,100) - #warp_m['.*warp_tanh.psi*'][:,0:1].constrain_fixed(1) - - #print(warp_m.checkgrad()) - #warp_m.plot() - #plt.show() + warp_m['.*noise.variance.*'].constrain_fixed(0.1) - warp_m.optimize_restarts(parallel=True, robust=True) - #print(warp_m.checkgrad()) + m = GPy.models.GPRegression(X[:, None], Y[:, None]) + m.optimize_restarts(parallel=False, robust=True, num_restarts=5) + warp_m.optimize_restarts(parallel=False, robust=True, num_restarts=5) print(warp_m) print(warp_m['.*warp.*']) warp_m.predict_in_warped_space = False warp_m.plot() warp_m.predict_in_warped_space = True warp_m.plot() + m.plot() warp_f.plot(X.min()-10, X.max()+10) + + + + plt.show() diff --git a/GPy/util/warping_functions.py b/GPy/util/warping_functions.py index 1617283c..516103bf 100644 --- a/GPy/util/warping_functions.py +++ b/GPy/util/warping_functions.py @@ -196,7 +196,11 @@ class TanhWarpingFunction_d(WarpingFunction): z = z.copy() if y is None: - y = np.ones_like(z) + # The idea here is to initialize y with +1 where + # z is positive and -1 where it is negative. + # For negative z, Newton-Raphson diverges + # if we initialize y with a positive value (and vice-versa). + y = ((z > 0) * 1.) - (z <= 0) it = 0 update = np.inf