(wpgs) fixing newton-raphson for f_inv and fixing plotting stuff

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):
+                median=False, deg_gauss_hermite=100, likelihood=None):
-        # normalize X values
+        """
-        # Xnew = (Xnew.copy() - self._Xoffset) / self._Xscale
+        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__':

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
+    return Xnew, xx, yy, xmin, xmax

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.sin(X) + np.random.normal(0,0.2,151)
-        Y = np.exp(Y) - 5
+        Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y])
-        #Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y]) + 0
         
-        #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['.*noise.variance.*'].constrain_fixed(0.1)
-        #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.optimize_restarts(parallel=True, robust=True)
+        m = GPy.models.GPRegression(X[:, None], Y[:, None])
-        #print(warp_m.checkgrad())
+        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()
 
         

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