diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py
index 487a324b..8ab00e10 100644
--- a/GPy/kern/constructors.py
+++ b/GPy/kern/constructors.py
@@ -296,5 +296,5 @@ def independent_outputs(k):
     """
     for sl in k.input_slices:
         assert (sl.start is None) and (sl.stop is None), "cannot adjust input slices! (TODO)"
-    parts = [independent_outputs.IndependentOutputs(p) for p in k.parts]
-    return kern(k.input_dim+1,parts)
+    new_parts = [parts.independent_outputs.IndependentOutputs(p) for p in k.parts]
+    return kern(k.input_dim+1,new_parts)
diff --git a/GPy/models/warped_gp.py b/GPy/models/warped_gp.py
index dc6eeb46..260139a6 100644
--- a/GPy/models/warped_gp.py
+++ b/GPy/models/warped_gp.py
@@ -19,7 +19,9 @@ class WarpedGP(GP):
             self.warping_function = TanhWarpingFunction_d(warping_terms)
             self.warping_params = (np.random.randn(self.warping_function.n_terms * 3 + 1,) * 1)
 
-        Y = self._scale_data(Y)
+        self.scale_data = False
+        if self.scale_data:
+            Y = self._scale_data(Y)
         self.has_uncertain_inputs = False
         self.Y_untransformed = Y.copy()
         self.predict_in_warped_space = False
@@ -87,11 +89,19 @@ class WarpedGP(GP):
     def plot_warping(self):
         self.warping_function.plot(self.warping_params, self.Y_untransformed.min(), self.Y_untransformed.max())
 
-    def _raw_predict(self, *args, **kwargs):
-        mu, var = GP._raw_predict(self, *args, **kwargs)
+    def predict(self, Xnew, which_parts='all', full_cov=False, pred_init=None):
+        # normalize X values
+        Xnew = (Xnew.copy() - self._Xoffset) / self._Xscale
+        mu, var = GP._raw_predict(self, Xnew, full_cov=full_cov, which_parts=which_parts)
+
+        # now push through likelihood
+        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov)
 
         if self.predict_in_warped_space:
-            mu = self.warping_function.f_inv(mu, self.warping_params)
+            mean = self.warping_function.f_inv(mean, self.warping_params, y=pred_init)
             var = self.warping_function.f_inv(var, self.warping_params)
-            mu = self._unscale_data(mu)
-        return mu, var
+
+        if self.scale_data:
+            mean = self._unscale_data(mean)
+        
+        return mean, var, _025pm, _975pm
diff --git a/GPy/util/warping_functions.py b/GPy/util/warping_functions.py
index 98d4d2b7..f36805a9 100644
--- a/GPy/util/warping_functions.py
+++ b/GPy/util/warping_functions.py
@@ -185,7 +185,7 @@ class TanhWarpingFunction_d(WarpingFunction):
         return z
 
 
-    def f_inv(self, z, psi, max_iterations = 1000):
+    def f_inv(self, z, psi, max_iterations=1000, y=None):
         """
         calculate the numerical inverse of f
 
@@ -195,7 +195,9 @@ class TanhWarpingFunction_d(WarpingFunction):
         """
 
         z = z.copy()
-        y = np.ones_like(z)
+        if y is None:
+            y = np.ones_like(z)
+            
         it = 0
         update = np.inf
 
@@ -205,7 +207,7 @@ class TanhWarpingFunction_d(WarpingFunction):
             it += 1
         if it == max_iterations:
             print "WARNING!!! Maximum number of iterations reached in f_inv "
-            
+
         return y