diff --git a/GPy/inference/latent_function_inference/posterior.py b/GPy/inference/latent_function_inference/posterior.py
index 56c0f388..8bfc1179 100644
--- a/GPy/inference/latent_function_inference/posterior.py
+++ b/GPy/inference/latent_function_inference/posterior.py
@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-from ...util.linalg import pdinv, dpotrs, tdot
+from ...util.linalg import pdinv, dpotrs, tdot, dtrtrs
 
 class Posterior(object):
     """
@@ -26,19 +26,19 @@ class Posterior(object):
         cov : the posterior covariance
 
         Not all of the above need to be supplied! You *must* supply:
-          
+
           log_marginal
           dL_dK
           dL_dtheta_lik
           K (for lazy computation)
 
        You may supply either:
-         
+         cc
           woodbury_chol
           woodbury_vector
 
         Or:
-          
+
           mean
           cov
           K_chol (for lazy computation)
@@ -46,7 +46,6 @@ class Posterior(object):
         Of course, you can supply more than that, but this class will lazily compute all other quantites on demand.
 
         From the supplied quantities, all of the others will be computed on demand (lazy computation)
-          
         """
         #obligatory
         self.log_marginal = log_marginal
@@ -80,7 +79,7 @@ class Posterior(object):
     @property
     def covariance(self):
         if self._covariance is None:
-            LiK, _ = dpotrs(self._woodbury_chol, self._K)
+            LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
             self._covariance = self._K - tdot(LiK.T)
         return self._covariance
 
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 988301f6..e3f21f14 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -80,6 +80,13 @@ class Gaussian(Likelihood):
         Z_hat = 1./np.sqrt(2.*np.pi*sum_var)*np.exp(-.5*(data_i - v_i/tau_i)**2./sum_var)
         return Z_hat, mu_hat, sigma2_hat
 
+    def predictive_values(self, mu, var, full_cov=False):
+        if full_cov:
+            low, up = mu - np.diag(var)[:,None], mu + np.diag(var)[:,None]
+        else:
+            low, up = mu - var, mu + var
+        return mu, var, low, up
+
     def predictive_mean(self, mu, sigma):
         #new_sigma2 = self.predictive_variance(mu, sigma)
         #return new_sigma2*(mu/sigma**2 + self.gp_link.transf(mu)/self.variance)