diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index 445f9ea8..b3f7115a 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -190,7 +190,7 @@ class GPBase(Model):
                 upper = m + 2*np.sqrt(v)
                 Y = self.Y
             else:
-                m, v, lower, upper = self.predict(Xgrid, which_parts=which_parts) #Compute the exact mean
+                m, v, lower, upper = self.predict(Xgrid, which_parts=which_parts)
                 Y = self.Y
             for d in which_data_ycols:
                 gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol)
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index e3f21f14..8f2d2225 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -54,7 +54,7 @@ class Gaussian(Likelihood):
 
     def _gradients(self, partial):
         """
-        Return the derivative of the log marginal likelihood wrt self.variance, 
+        Return the derivative of the log marginal likelihood wrt self.variance,
         given the appropriate partial derivative
         """
         return np.sum(partial)
@@ -82,9 +82,13 @@ class Gaussian(Likelihood):
 
     def predictive_values(self, mu, var, full_cov=False):
         if full_cov:
-            low, up = mu - np.diag(var)[:,None], mu + np.diag(var)[:,None]
+            var += np.eye(var.shape[0])*self.variance
+            d = 2*np.sqrt(np.diag(var))
+            low, up = mu - d, mu + d
         else:
-            low, up = mu - var, mu + var
+            var += self.variance
+            d = 2*np.sqrt(var)
+            low, up = mu - d, mu + d
         return mu, var, low, up
 
     def predictive_mean(self, mu, sigma):