[beckdaniel] merge in warped gp changes

GPy/models/warped_gp.py

@@ -16,16 +16,16 @@ class WarpedGP(GP):
 
         if warping_function == None:
             self.warping_function = TanhWarpingFunction_d(warping_terms)
-            self.warping_params = (np.random.randn(self.warping_function.n_terms * 3 + 1,) * 1)
+            self.warping_params = (np.random.randn(self.warping_function.n_terms * 3 + 1) * 1)
         else:
             self.warping_function = warping_function
 
         self.scale_data = False
         if self.scale_data:
             Y = self._scale_data(Y)
-        self.has_uncertain_inputs = False
+        #self.has_uncertain_inputs = False
         self.Y_untransformed = Y.copy()
-        self.predict_in_warped_space = False
+        self.predict_in_warped_space = True
         likelihood = likelihoods.Gaussian()
 
         GP.__init__(self, X, self.transform_data(), likelihood=likelihood, kernel=kernel)
@@ -56,7 +56,6 @@ class WarpedGP(GP):
         self.warping_function.psi.gradient[:] = warping_grads[:, :-1]
         self.warping_function.d.gradient[:] = warping_grads[0, -1]
 
-
     def transform_data(self):
         Y = self.warping_function.f(self.Y_untransformed.copy()).copy()
         return Y
@@ -69,7 +68,35 @@ class WarpedGP(GP):
     def plot_warping(self):
         self.warping_function.plot(self.Y_untransformed.min(), self.Y_untransformed.max())
 
-    def predict(self, Xnew, which_parts='all', pred_init=None):
+    def _get_warped_term(self, mean, std, gh_samples, pred_init=None):
+        arg1 = gh_samples.dot(std.T) * np.sqrt(2)
+        arg2 = np.ones(shape=gh_samples.shape).dot(mean.T)
+        return self.warping_function.f_inv(arg1 + arg2, y=pred_init)
+
+    def _get_warped_mean(self, mean, std, pred_init=None, deg_gauss_hermite=100):
+        """
+        Calculate the warped mean by using Gauss-Hermite quadrature.
+        """
+        gh_samples, gh_weights = np.polynomial.hermite.hermgauss(deg_gauss_hermite)
+        gh_samples = gh_samples[:,None]
+        gh_weights = gh_weights[None,:]
+        return gh_weights.dot(self._get_warped_term(mean, std, gh_samples)) / np.sqrt(np.pi)
+
+    def _get_warped_variance(self, mean, std, pred_init=None, deg_gauss_hermite=100):
+        """
+        Calculate the warped variance by using Gauss-Hermite quadrature.
+        """
+        gh_samples, gh_weights = np.polynomial.hermite.hermgauss(deg_gauss_hermite)
+        gh_samples = gh_samples[:,None]
+        gh_weights = gh_weights[None,:]
+        arg1 = gh_weights.dot(self._get_warped_term(mean, std, gh_samples, 
+                                                    pred_init=pred_init) ** 2) / np.sqrt(np.pi)
+        arg2 = self._get_warped_mean(mean, std, pred_init=pred_init,
+                                     deg_gauss_hermite=deg_gauss_hermite)
+        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
         mu, var = GP._raw_predict(self, Xnew)
@@ -78,13 +105,48 @@ class WarpedGP(GP):
         mean, var = self.likelihood.predictive_values(mu, var)
 
         if self.predict_in_warped_space:
-            mean = self.warping_function.f_inv(mean,  y=pred_init)
-            var = self.warping_function.f_inv(var)
+            std = np.sqrt(var)
+            if median:
+                wmean = self.warping_function.f_inv(mean, y=pred_init)
+            else:
+                wmean = self._get_warped_mean(mean, std, pred_init=pred_init,
+                                              deg_gauss_hermite=deg_gauss_hermite).T
+            wvar = self._get_warped_variance(mean, std, pred_init=pred_init,
+                                             deg_gauss_hermite=deg_gauss_hermite).T
+        else:
+            wmean = mean
+            wvar = var
 
         if self.scale_data:
-            mean = self._unscale_data(mean)
+            pred = self._unscale_data(pred)
+
+        return wmean, wvar
+
+    def predict_quantiles(self, X, quantiles=(2.5, 97.5), Y_metadata=None):
+        """
+        Get the predictive quantiles around the prediction at X
+
+        :param X: The points at which to make a prediction
+        :type X: np.ndarray (Xnew x self.input_dim)
+        :param quantiles: tuple of quantiles, default is (2.5, 97.5) which is the 95% interval
+        :type quantiles: tuple
+        :returns: list of quantiles for each X and predictive quantiles for interval combination
+        :rtype: [np.ndarray (Xnew x self.input_dim), np.ndarray (Xnew x self.input_dim)]
+        """
+        m, v = self._raw_predict(X,  full_cov=False)
+        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)
 
-        return mean, var
 
 if __name__ == '__main__':
     X = np.random.randn(100, 1)