diff --git a/GPy/models/warped_gp.py b/GPy/models/warped_gp.py
index 06bc85cb..df947d3e 100644
--- a/GPy/models/warped_gp.py
+++ b/GPy/models/warped_gp.py
@@ -137,6 +137,7 @@ 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
diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py
index 0986c20e..8fe0b2c2 100644
--- a/GPy/plotting/matplot_dep/models_plots.py
+++ b/GPy/plotting/matplot_dep/models_plots.py
@@ -46,6 +46,8 @@ def plot_data(model, which_data_rows='all',
     #data
     X = model.X
     Y = model.Y
+    if isinstance(model, WarpedGP) and model.predict_in_warped_space:
+        Y = model.Y_untransformed
 
     #work out what the inputs are for plotting (1D or 2D)
     if visible_dims is None:
@@ -184,11 +186,12 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
             if isinstance(model, WarpedGP):
                 m, v = model.predict(Xgrid, full_cov=False, median=True, Y_metadata=Y_metadata, **predict_kw)
+                lower, upper = model.predict_quantiles(Xgrid)
                 #print np.concatenate((Xgrid, m), axis=1)
             else:
                 m, v = model.predict(Xgrid, full_cov=False, Y_metadata=Y_metadata, **predict_kw)
-            fmu, fv = model._raw_predict(Xgrid, full_cov=False, **predict_kw)
-            lower, upper = model.likelihood.predictive_quantiles(fmu, fv, (2.5, 97.5), Y_metadata=Y_metadata)
+                fmu, fv = model._raw_predict(Xgrid, full_cov=False, **predict_kw)
+                lower, upper = model.likelihood.predictive_quantiles(fmu, fv, (2.5, 97.5), Y_metadata=Y_metadata)