diff --git a/GPy/examples/ep_fix.py b/GPy/examples/ep_fix.py
index 3cb35663..d1747025 100644
--- a/GPy/examples/ep_fix.py
+++ b/GPy/examples/ep_fix.py
@@ -4,7 +4,7 @@
 
 """
 Simple Gaussian Processes classification 1D
-Probit likelihood
+probit likelihood
 """
 import pylab as pb
 import numpy as np
@@ -26,7 +26,7 @@ Y = np.hstack([np.ones(N/2),np.repeat(-1,N/2)])[:,None]
 kernel = GPy.kern.rbf(1)
 
 # Define likelihood
-distribution = GPy.likelihoods.likelihood_functions.Probit()
+distribution = GPy.likelihoods.likelihood_functions.probit()
 likelihood_object = GPy.likelihoods.EP(Y,distribution)
 
 # Model definition
diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py
index 420b138a..a88059b1 100644
--- a/GPy/likelihoods/EP.py
+++ b/GPy/likelihoods/EP.py
@@ -27,7 +27,7 @@ class EP(likelihood):
 
         #initial values for the GP variables
         self.Y = np.zeros((self.N,1))
-        self.variance = np.zeros((self.N,self.N))#np.eye(self.N)
+        self.covariance_matrix = np.eye(self.N)
         self.Z = 0
         self.YYT = None
 
@@ -50,8 +50,9 @@ class EP(likelihood):
         mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2
         self.Z = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant, aka Z_ep
 
-        self.Y =  mu_tilde[:,None]
-        self.precsion = self.tau_tilde[:,None]
+        self.Y = mu_tilde[:,None]
+        self.YYT = np.dot(self.Y,self.Y.T)
+        self.precision = self.tau_tilde
         self.covariance_matrix = np.diag(1./self.precision)
 
     def fit_full(self,K):
@@ -61,9 +62,11 @@ class EP(likelihood):
         """
         #Prior distribution parameters: p(f|X) = N(f|0,K)
 
+        self.tau_tilde = np.zeros(self.N)
+        self.v_tilde = np.zeros(self.N)
         #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
         self.mu = np.zeros(self.N)
-        self.Sigma = K.copy() - self.variance.copy()
+        self.Sigma = K.copy()
 
         """
         Initial values - Cavity distribution parameters:
diff --git a/GPy/models/GP.py b/GPy/models/GP.py
index ae192618..e64da2c9 100644
--- a/GPy/models/GP.py
+++ b/GPy/models/GP.py
@@ -65,7 +65,7 @@ class GP(model):
         self.likelihood._set_params(p[self.kern.Nparam:])
 
         self.K = self.kern.K(self.X,slices1=self.Xslices)
-        self.K += self.likelihood.variance
+        self.K += self.likelihood.covariance_matrix
 
         self.Ki, self.L, self.Li, self.K_logdet = pdinv(self.K)
 
@@ -90,7 +90,7 @@ class GP(model):
         For a Gaussian (or direct: TODO) likelihood, no iteration is required:
         this function does nothing
         """
-        self.likelihood.fit_full(self.kern.compute(self.X))
+        self.likelihood.fit_full(self.kern.K(self.X))
         self._set_params(self._get_params()) # update the GP
 
     def _model_fit_term(self):