diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py
index f49ed275..9c55e5f7 100644
--- a/GPy/likelihoods/EP.py
+++ b/GPy/likelihoods/EP.py
@@ -33,7 +33,9 @@ class EP(likelihood):
         self.Z = 0
         self.YYT = None
 
-    def predictive_values(self,mu,var):
+    def predictive_values(self,mu,var,full_cov):
+        if full_cov:
+            raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood"
         return self.likelihood_function.predictive_values(mu,var)
 
     def _get_params(self):
diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py
index fc36a63a..25d12491 100644
--- a/GPy/likelihoods/Gaussian.py
+++ b/GPy/likelihoods/Gaussian.py
@@ -43,14 +43,24 @@ class Gaussian(likelihood):
         self.covariance_matrix = np.eye(self.N)*self._variance
         self.precision = 1./self._variance
 
-    def predictive_values(self,mu,var):
+    def predictive_values(self,mu,var, full_cov):
         """
         Un-normalize the prediction and add the likelihood variance, then return the 5%, 95% interval
         """
         mean = mu*self._std + self._mean
-        true_var = (var + self._variance)*self._std**2
-        _5pc = mean + - 2.*np.sqrt(true_var)
-        _95pc = mean + 2.*np.sqrt(true_var)
+        if full_cov:
+            if self.D >1:
+                raise NotImplementedError, "TODO"
+                #Note. for D>1, we need to re-normalise all the outputs independently.
+                # This will mess up computations of diag(true_var), below.
+                #note that the upper, lower quantiles should be the same shape as mean
+            true_var = (var + np.eye(var.shape[0])*self._variance)*self._std**2
+            _5pc = mean + - 2.*np.sqrt(np.diag(true_var))
+            _95pc = mean + 2.*np.sqrt(np.diag(true_var))
+        else:
+            true_var = (var + self._variance)*self._std**2
+            _5pc = mean + - 2.*np.sqrt(true_var)
+            _95pc = mean + 2.*np.sqrt(true_var)
         return mean, true_var, _5pc, _95pc
 
     def fit_full(self):
diff --git a/GPy/models/GP.py b/GPy/models/GP.py
index 238e7cd7..53ba1183 100644
--- a/GPy/models/GP.py
+++ b/GPy/models/GP.py
@@ -140,7 +140,7 @@ class GP(model):
         KiKx = np.dot(self.Ki,Kx)
         if full_cov:
             Kxx = self.kern.K(_Xnew, slices1=slices,slices2=slices)
-            var = Kxx - np.dot(KiKx.T,Kx) #NOTE this won't work for plotting
+            var = Kxx - np.dot(KiKx.T,Kx)
         else:
             Kxx = self.kern.Kdiag(_Xnew, slices=slices)
             var = Kxx - np.sum(np.multiply(KiKx,Kx),0)
@@ -179,7 +179,7 @@ class GP(model):
         mu, var = self._raw_predict(Xnew, slices, full_cov)
 
         #now push through likelihood TODO
-        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var)
+        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov)
 
         return mean, var, _025pm, _975pm