diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index b6e4ebc0..cb968520 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -173,7 +173,8 @@ class GPBase(Model):
                 upper = m + 2*np.sqrt(v)
                 Y = self.likelihood.Y
             else:
-                m, v, lower, upper = self.predict(Xgrid, which_parts=which_parts)
+                m, v, lower, upper = self.predict(Xgrid, which_parts=which_parts,sampling=False) #Compute the exact mean
+                m_, v_, lower, upper = self.predict(Xgrid, which_parts=which_parts,sampling=True,num_samples=15000) #Apporximate the percentiles
                 Y = self.likelihood.data
             for d in which_data_ycols:
                 gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol)
@@ -210,7 +211,7 @@ class GPBase(Model):
                 m, _ = self._raw_predict(Xgrid, which_parts=which_parts)
                 Y = self.likelihood.Y
             else:
-                m, _, _, _ = self.predict(Xgrid, which_parts=which_parts,num_samples=100) #FIXME we need a balance between accuracy and speed to define num_samples
+                m, _, _, _ = self.predict(Xgrid, which_parts=which_parts,sampling=False)
                 Y = self.likelihood.data
             for d in which_data_ycols:
                 m_d = m[:,d].reshape(resolution, resolution).T
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 85c028b4..c12d8e6d 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -69,7 +69,7 @@ class Gaussian(likelihood):
             self.covariance_matrix = np.eye(self.N) * x
             self._variance = x
 
-    def predictive_values(self, mu, var, full_cov):
+    def predictive_values(self, mu, var, full_cov, **likelihood_args):
         """
         Un-normalize the prediction and add the likelihood variance, then return the 5%, 95% interval
         """