diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py
index a5084cc0..fc36a63a 100644
--- a/GPy/likelihoods/Gaussian.py
+++ b/GPy/likelihoods/Gaussian.py
@@ -51,7 +51,7 @@ class Gaussian(likelihood):
         true_var = (var + self._variance)*self._std**2
         _5pc = mean + - 2.*np.sqrt(true_var)
         _95pc = mean + 2.*np.sqrt(true_var)
-        return mean, _5pc, _95pc
+        return mean, true_var, _5pc, _95pc
 
     def fit_full(self):
         """
diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py
index 6ec57c07..c1d9585e 100644
--- a/GPy/likelihoods/likelihood.py
+++ b/GPy/likelihoods/likelihood.py
@@ -25,11 +25,16 @@ class likelihood:
     def _get_param_names(self):
         raise NotImplementedError
 
-    def _set_params(self,x):
+    def _set_params(self, x):
         raise NotImplementedError
 
     def fit(self):
         raise NotImplementedError
 
-    def _gradients(self,partial):
+    def _gradients(self, partial):
         raise NotImplementedError
+
+    def predictive_values(self, mu, var):
+        raise NotImplementedError
+
+
diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py
index 3e2a0361..4b8e7013 100644
--- a/GPy/likelihoods/likelihood_functions.py
+++ b/GPy/likelihoods/likelihood_functions.py
@@ -48,14 +48,14 @@ class probit(likelihood_function):
 
     def predictive_values(self,mu,var):
         """
-        Compute  mean, and conficence interval (percentiles 5 and 95) of the  prediction
+        Compute  mean, variance and conficence interval (percentiles 5 and 95) of the  prediction
         """
         mu = mu.flatten()
         var = var.flatten()
         mean = stats.norm.cdf(mu/np.sqrt(1+var))
         p_025 = np.zeros(mu.shape)
         p_975 = np.ones(mu.shape)
-        return mean, p_025, p_975
+        return mean, np.nan*var, p_025, p_975 # TODO: better values here (mean is okay)
 
 class Poisson(likelihood_function):
     """
@@ -131,4 +131,4 @@ class Poisson(likelihood_function):
         tmp = stats.poisson.ppf(np.array([.025,.975]),mean)
         p_025 = tmp[:,0]
         p_975 = tmp[:,1]
-        return mean,p_025,p_975
+        return mean,np.nan*mean,p_025,p_975 # better variance here TODO
diff --git a/GPy/models/GP.py b/GPy/models/GP.py
index 796ab7d6..238e7cd7 100644
--- a/GPy/models/GP.py
+++ b/GPy/models/GP.py
@@ -179,7 +179,7 @@ class GP(model):
         mu, var = self._raw_predict(Xnew, slices, full_cov)
 
         #now push through likelihood TODO
-        mean, _025pm, _975pm = self.likelihood.predictive_values(mu, var)
+        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var)
 
         return mean, var, _025pm, _975pm