diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index 21b435e7..452167ce 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -83,7 +83,7 @@ def coregionalisation_toy2(optim_iters=100):
     Y = np.vstack((Y1,Y2))
 
     k1 = GPy.kern.rbf(1) + GPy.kern.bias(1)
-    k2 = GPy.kern.Coregionalise(2,1)
+    k2 = GPy.kern.coregionalise(2,1)
     k = k1.prod(k2,tensor=True)
     m = GPy.models.GPRegression(X,Y,kernel=k)
     m.constrain_fixed('.*rbf_var',1.)
@@ -114,7 +114,7 @@ def coregionalisation_toy(optim_iters=100):
     Y = np.vstack((Y1,Y2))
 
     k1 = GPy.kern.rbf(1)
-    k2 = GPy.kern.Coregionalise(2,2)
+    k2 = GPy.kern.coregionalise(2,2)
     k = k1.prod(k2,tensor=True)
     m = GPy.models.GPRegression(X,Y,kernel=k)
     m.constrain_fixed('.*rbf_var',1.)
@@ -149,7 +149,7 @@ def coregionalisation_sparse(optim_iters=100):
     Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None]))
 
     k1 = GPy.kern.rbf(1)
-    k2 = GPy.kern.Coregionalise(2,2)
+    k2 = GPy.kern.coregionalise(2,2)
     k = k1.prod(k2,tensor=True) + GPy.kern.white(2,0.001)
 
     m = GPy.models.SparseGPRegression(X,Y,kernel=k,Z=Z)
diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py
index e058de79..4932bd40 100644
--- a/GPy/likelihoods/__init__.py
+++ b/GPy/likelihoods/__init__.py
@@ -1,6 +1,5 @@
 from ep import EP
 from gaussian import Gaussian
 # TODO: from Laplace import Laplace
-import likelihood_functions as functions
-import binomial_likelihood
-import poisson_likelihood
+from constructors import *
+
diff --git a/GPy/likelihoods/binomial_likelihood.py b/GPy/likelihoods/binomial_likelihood.py
deleted file mode 100644
index 15b8067e..00000000
--- a/GPy/likelihoods/binomial_likelihood.py
+++ /dev/null
@@ -1,123 +0,0 @@
-# Copyright (c) 2012, 2013 Ricardo Andrade
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-from scipy import stats,special
-import scipy as sp
-import pylab as pb
-from ..util.plot import gpplot
-from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import link_functions
-from likelihood_functions import LikelihoodFunction
-
-class Binomial(LikelihoodFunction):
-    """
-    Probit likelihood
-    Y is expected to take values in {-1,1}
-    -----
-    $$
-    L(x) = \\Phi (Y_i*f_i)
-    $$
-    """
-    def __init__(self,link=None):
-        self.discrete = True
-        self.support_limits = (0,1)
-
-        if not link:
-            link = link_functions.Probit
-        if isinstance(link,link_functions.Probit):
-            self.analytical_moments = True
-        else:
-            self.analytical_moments = False
-        super(Binomial, self).__init__(link)
-
-    def _preprocess_values(self,Y):
-        """
-        Check if the values of the observations correspond to the values
-        assumed by the likelihood function.
-
-        ..Note:: Binary classification algorithm works better with classes {-1,1}
-        """
-        Y_prep = Y.copy()
-        Y1 = Y[Y.flatten()==1].size
-        Y2 = Y[Y.flatten()==0].size
-        assert Y1 + Y2 == Y.size, 'Binomial likelihood is meant to be used only with outputs in {0,1}.'
-        Y_prep[Y.flatten() == 0] = -1
-        return Y_prep
-
-    def _moments_match_analytical(self,data_i,tau_i,v_i):
-        """
-        Moments match of the marginal approximation in EP algorithm
-
-        :param i: number of observation (int)
-        :param tau_i: precision of the cavity distribution (float)
-        :param v_i: mean/variance of the cavity distribution (float)
-        """
-        z = data_i*v_i/np.sqrt(tau_i**2 + tau_i)
-        Z_hat = std_norm_cdf(z)
-        phi = std_norm_pdf(z)
-        mu_hat = v_i/tau_i + data_i*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
-        sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
-        return Z_hat, mu_hat, sigma2_hat
-
-    def _predictive_mean_analytical(self,mu,sigma):
-        return stats.norm.cdf(mu/np.sqrt(1+sigma**2))
-
-    def _mass(self,gp,obs):
-        #NOTE obs must be in {0,1}
-        p = self.link.inv_transf(gp)
-        return p**obs * (1.-p)**(1.-obs)
-
-    def _nlog_mass(self,gp,obs):
-        p = self.link.inv_transf(gp)
-        return obs*np.log(p) + (1.-obs)*np.log(1-p)
-
-    def _dnlog_mass_dgp(self,gp,obs):
-        p = self.link.inv_transf(gp)
-        dp = self.link.dinv_transf_df(gp)
-        return obs/p * dp - (1.-obs)/(1.-p) * dp
-
-    def _d2nlog_mass_dgp2(self,gp,obs):
-        p = self.link.inv_transf(gp)
-        return (obs/p + (1.-obs)/(1.-p))*self.lind.d2inv_transf_df(gp) + ((1.-obs)/(1.-p)**2-obs/p**2)*self.link.dinv_transf_df(gp)
-
-    def _mean(self,gp):
-        """
-        Mass (or density) function
-        """
-        return self.link.inv_transf(gp)
-
-    def _dmean_dgp(self,gp):
-        return self.link.dinv_transf_df(gp)
-
-    def _d2mean_dgp2(self,gp):
-        return self.link.d2inv_transf_df2(gp)
-
-    def _variance(self,gp):
-        """
-        Mass (or density) function
-        """
-        p = self.link.inv_transf(gp)
-        return p*(1-p)
-
-    def _dvariance_dgp(self,gp):
-        return self.link.dinv_transf_df(gp)*(1. - 2.*self.link.inv_transf(gp))
-
-    def _d2variance_dgp2(self,gp):
-        return self.link.d2inv_transf_df2(gp)*(1. - 2.*self.link.inv_transf(gp)) - 2*self.link.dinv_transf_df(gp)**2
-
-    """
-    def predictive_values(self,mu,var): #TODO remove
-        mu = mu.flatten()
-        var = var.flatten()
-        #mean = stats.norm.cdf(mu/np.sqrt(1+var))
-        mean = self._predictive_mean_analytical(mu,np.sqrt(var))
-        norm_025 = [stats.norm.ppf(.025,m,v) for m,v in zip(mu,var)]
-        norm_975 = [stats.norm.ppf(.975,m,v) for m,v in zip(mu,var)]
-        #p_025 = stats.norm.cdf(norm_025/np.sqrt(1+var))
-        #p_975 = stats.norm.cdf(norm_975/np.sqrt(1+var))
-        p_025 = self._predictive_mean_analytical(norm_025,np.sqrt(var))
-        p_975 = self._predictive_mean_analytical(norm_975,np.sqrt(var))
-        return mean[:,None], np.nan*var, p_025[:,None], p_975[:,None] # TODO: var
-    """
diff --git a/GPy/likelihoods/constructors.py b/GPy/likelihoods/constructors.py
new file mode 100644
index 00000000..0b995894
--- /dev/null
+++ b/GPy/likelihoods/constructors.py
@@ -0,0 +1,42 @@
+# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from likelihood_functions import LikelihoodFunction
+import noise_models
+import link_functions
+
+def binomial(link=None):
+    """
+    Construct a binomial likelihood
+
+    :param link: a GPy link function
+    """
+    #self.discrete = True
+    #self.support_limits = (0,1)
+
+    if link is None:
+        link = link_functions.Probit()
+    else:
+        assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
+
+    if isinstance(link,link_functions.Probit):
+        analytical_moments = True
+    else:
+        analytical_moments = False
+    return noise_models.binomial_likelihood.Binomial(link,analytical_moments)
+
+
+def poisson(link=None):
+    """
+    Construct a Poisson likelihood
+
+    :param link: a GPy link function
+    """
+    if link is None:
+        link = link_functions.Log_ex_1()
+    else:
+        assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
+    #assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
+    analytical_moments = False
+    return noise_models.poisson_likelihood.Poisson(link,analytical_moments)
diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py
deleted file mode 100644
index 9b132d72..00000000
--- a/GPy/likelihoods/likelihood_functions.py
+++ /dev/null
@@ -1,348 +0,0 @@
-# Copyright (c) 2012, 2013 Ricardo Andrade
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-from scipy import stats,special
-import scipy as sp
-import pylab as pb
-from ..util.plot import gpplot
-from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import link_functions
-
-
-class LikelihoodFunction(object):
-    """
-    Likelihood class for doing Expectation propagation
-
-    :param Y: observed output (Nx1 numpy.darray)
-    ..Note:: Y values allowed depend on the LikelihoodFunction used
-    """
-    def __init__(self,link):
-        #assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME
-        self.link = link
-
-        if self.analytical_moments:
-            self.moments_match = self._moments_match_analytical
-            self.predictive_mean = self._predictive_mean_analytical
-        else:
-            self.moments_match = self._moments_match_numerical
-            self.predictive_mean = self._predictive_mean_numerical
-
-    def _preprocess_values(self,Y):
-        """
-        In case it is needed, this function assess the output values or makes any pertinent transformation on them.
-
-        :param Y: observed output (Nx1 numpy.darray)
-        """
-        return Y
-
-    def _product(self,gp,obs,mu,sigma):
-        """
-        Product between the cavity distribution and a likelihood factor.
-
-        :param gp: latent variable
-        :param obs: observed output
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return stats.norm.pdf(gp,loc=mu,scale=sigma) * self._mass(gp,obs)
-
-    def _nlog_product_scaled(self,gp,obs,mu,sigma):
-        """
-        Negative log-product between the cavity distribution and a likelihood factor.
-        ..Note:: The constant term in the Gaussian distribution is ignored.
-
-        :param gp: latent variable
-        :param obs: observed output
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return .5*((gp-mu)/sigma)**2 + self._nlog_mass(gp,obs)
-
-    def _dnlog_product_dgp(self,gp,obs,mu,sigma):
-        """
-        Derivative wrt latent variable of the log-product between the cavity distribution and a likelihood factor.
-
-        :param gp: latent variable
-        :param obs: observed output
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return (gp - mu)/sigma**2 + self._dnlog_mass_dgp(gp,obs)
-
-    def _d2nlog_product_dgp2(self,gp,obs,mu,sigma):
-        """
-        Second derivative wrt latent variable of the log-product between the cavity distribution and a likelihood factor.
-
-        :param gp: latent variable
-        :param obs: observed output
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return 1./sigma**2 + self._d2nlog_mass_dgp2(gp,obs)
-
-    def _product_mode(self,obs,mu,sigma):
-        """
-        Newton's CG method to find the mode in _product (cavity x likelihood factor).
-
-        :param obs: observed output
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return sp.optimize.fmin_ncg(self._nlog_product_scaled,x0=mu,fprime=self._dnlog_product_dgp,fhess=self._d2nlog_product_dgp2,args=(obs,mu,sigma))
-
-    def _moments_match_analytical(self,obs,tau,v):
-        """
-        If available, this function computes the moments analytically.
-        """
-        pass
-
-    def _moments_match_numerical(self,obs,tau,v):
-        """
-        Lapace approximation to calculate the moments.
-
-        :param obs: observed output
-        :param tau: cavity distribution 1st natural parameter (precision)
-        :param v: cavity distribution 2nd natural paramenter (mu*precision)
-        """
-        mu = v/tau
-        mu_hat = self._product_mode(obs,mu,np.sqrt(1./tau))
-        sigma2_hat = 1./(tau + self._d2nlog_mass_dgp2(mu_hat,obs))
-        Z_hat = np.exp(-.5*tau*(mu_hat-mu)**2) * self._mass(mu_hat,obs)*np.sqrt(tau*sigma2_hat)
-        return Z_hat,mu_hat,sigma2_hat
-
-    def _nlog_conditional_mean_scaled(self,gp,mu,sigma):
-        """
-        Negative logarithm of the l.v.'s predictive distribution times the output's mean given the l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-
-        ..Note:: This function helps computing E(Y_star) = E(E(Y_star|f_star))
-        """
-        return .5*((gp - mu)/sigma)**2 - np.log(self._mean(gp))
-
-    def _dnlog_conditional_mean_dgp(self,gp,mu,sigma):
-        """
-        Derivative of _nlog_conditional_mean_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return (gp - mu)/sigma**2 - self._dmean_dgp(gp)/self._mean(gp)
-
-    def _d2nlog_conditional_mean_dgp2(self,gp,mu,sigma):
-        """
-        Second derivative of _nlog_conditional_mean_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return 1./sigma**2 - self._d2mean_dgp2(gp)/self._mean(gp) + (self._dmean_dgp(gp)/self._mean(gp))**2
-
-    def _nlog_exp_conditional_variance_scaled(self,gp,mu,sigma):
-        """
-        Negative logarithm of the l.v.'s predictive distribution times the output's variance given the l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-
-        ..Note:: This function helps computing E(V(Y_star|f_star))
-        """
-        return .5*((gp - mu)/sigma)**2 - np.log(self._variance(gp))
-
-    def _dnlog_exp_conditional_variance_dgp(self,gp,mu,sigma):
-        """
-        Derivative of _nlog_exp_conditional_variance_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return (gp - mu)/sigma**2 - self._dvariance_dgp(gp)/self._variance(gp)
-
-    def _d2nlog_exp_conditional_variance_dgp2(self,gp,mu,sigma):
-        """
-        Second derivative of _nlog_exp_conditional_variance_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return 1./sigma**2 - self._d2variance_dgp2(gp)/self._variance(gp) + (self._dvariance_dgp(gp)/self._variance(gp))**2
-
-    def _nlog_exp_conditional_mean_sq_scaled(self,gp,mu,sigma):
-        """
-        Negative logarithm of the l.v.'s predictive distribution times the output's mean squared given the l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-
-        ..Note:: This function helps computing E( E(Y_star|f_star)**2 )
-        """
-        return .5*((gp - mu)/sigma)**2 - 2*np.log(self._mean(gp))
-
-    def _dnlog_exp_conditional_mean_sq_dgp(self,gp,mu,sigma):
-        """
-        Derivative of _nlog_exp_conditional_mean_sq_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return (gp - mu)/sigma**2 - 2*self._dmean_dgp(gp)/self._mean(gp)
-
-    def _d2nlog_exp_conditional_mean_sq_dgp2(self,gp,mu,sigma):
-        """
-        Second derivative of _nlog_exp_conditional_mean_sq_scaled wrt. l.v.
-
-        :param gp: latent variable
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        return 1./sigma**2 - 2*( self._d2mean_dgp2(gp)/self._mean(gp) - (self._dmean_dgp(gp)/self._mean(gp))**2 )
-
-    def _predictive_mean_analytical(self,mu,sigma):
-        """
-        If available, this function computes the predictive mean analytically.
-        """
-        pass
-
-    def _predictive_mean_numerical(self,mu,sigma):
-        """
-        Laplace approximation to the predictive mean: E(Y_star) = E( E(Y_star|f_star) )
-
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        maximum = sp.optimize.fmin_ncg(self._nlog_conditional_mean_scaled,x0=self._mean(mu),fprime=self._dnlog_conditional_mean_dgp,fhess=self._d2nlog_conditional_mean_dgp2,args=(mu,sigma))
-        mean = np.exp(-self._nlog_conditional_mean_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_conditional_mean_dgp2(maximum,mu,sigma))*sigma)
-        """
-        pb.figure()
-        x = np.array([mu + step*sigma for step in np.linspace(-7,7,100)])
-        f = np.array([np.exp(-self._nlog_conditional_mean_scaled(xi,mu,sigma))/np.sqrt(2*np.pi*sigma**2) for xi in x])
-        pb.plot(x,f,'b-')
-        sigma2 = 1./self._d2nlog_conditional_mean_dgp2(maximum,mu,sigma)
-        f2 = np.exp(-.5*(x-maximum)**2/sigma2)/np.sqrt(2*np.pi*sigma2)
-        k = np.exp(-self._nlog_conditional_mean_scaled(maximum,mu,sigma))*np.sqrt(sigma2)/np.sqrt(sigma**2)
-        pb.plot(x,f2*mean,'r-')
-        pb.vlines(maximum,0,f.max())
-        """
-        return mean
-
-    def _predictive_mean_sq(self,mu,sigma):
-        """
-        Laplace approximation to the predictive mean squared: E(Y_star**2) = E( E(Y_star|f_star)**2 )
-
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        """
-        maximum = sp.optimize.fmin_ncg(self._nlog_exp_conditional_mean_sq_scaled,x0=self._mean(mu),fprime=self._dnlog_exp_conditional_mean_sq_dgp,fhess=self._d2nlog_exp_conditional_mean_sq_dgp2,args=(mu,sigma))
-        mean_squared = np.exp(-self._nlog_exp_conditional_mean_sq_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_mean_sq_dgp2(maximum,mu,sigma))*sigma)
-        return mean_squared
-
-    def predictive_variance(self,mu,sigma,predictive_mean=None):
-        """
-        Laplace approximation to the predictive variance: V(Y_star) = E( V(Y_star|f_star) ) + V( E(Y_star|f_star) )
-
-        :param mu: cavity distribution mean
-        :param sigma: cavity distribution standard deviation
-        :predictive_mean: output's predictive mean, if None _predictive_mean function will be called.
-        """
-        # E( V(Y_star|f_star) )
-        maximum = sp.optimize.fmin_ncg(self._nlog_exp_conditional_variance_scaled,x0=self._variance(mu),fprime=self._dnlog_exp_conditional_variance_dgp,fhess=self._d2nlog_exp_conditional_variance_dgp2,args=(mu,sigma))
-        exp_var = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_variance_dgp2(maximum,mu,sigma))*sigma)
-
-        """
-        pb.figure()
-        x = np.array([mu + step*sigma for step in np.linspace(-7,7,100)])
-        f = np.array([np.exp(-self._nlog_exp_conditional_variance_scaled(xi,mu,sigma))/np.sqrt(2*np.pi*sigma**2) for xi in x])
-        pb.plot(x,f,'b-')
-        sigma2 = 1./self._d2nlog_exp_conditional_variance_dgp2(maximum,mu,sigma)
-        f2 = np.exp(-.5*(x-maximum)**2/sigma2)/np.sqrt(2*np.pi*sigma2)
-        k = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))*np.sqrt(sigma2)/np.sqrt(sigma**2)
-        pb.plot(x,f2*exp_var,'r--')
-        pb.vlines(maximum,0,f.max())
-        """
-
-        #V( E(Y_star|f_star) ) =  E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star)**2 )
-        exp_exp2 = self._predictive_mean_sq(mu,sigma)
-        if predictive_mean is None:
-            predictive_mean = self.predictive_mean(mu,sigma)
-        var_exp = exp_exp2 - predictive_mean**2
-        return exp_var + var_exp
-
-    def _nlog_joint_predictive_scaled(self,x,mu,sigma):
-        """
-        Negative logarithm of the joint predictive distribution (latent variable and output).
-
-        :param x: tuple (latent variable,output)
-        :param mu: latent variable's predictive mean
-        :param sigma: latent variable's predictive standard deviation
-        """
-        return self._nlog_product_scaled(x[0],x[1],mu,sigma)
-
-    def _gradient_nlog_joint_predictive(self,x,mu,sigma):
-        """
-        Gradient of _nlog_joint_predictive_scaled.
-
-        :param x: tuple (latent variable,output)
-        :param mu: latent variable's predictive mean
-        :param sigma: latent variable's predictive standard deviation
-        ..Note: Only avilable when the output is continuous
-        """
-        assert not self.discrete, "Gradient not available for discrete outputs."
-        return np.array((self._dnlog_product_dgp(gp=x[0],obs=x[1],mu=mu,sigma=sigma),self._dnlog_mass_dobs(obs=x[1],gp=x[0])))
-
-    def _hessian_nlog_joint_predictive(self,x,mu,sigma):
-        """
-        Hessian of _nlog_joint_predictive_scaled.
-
-        :param x: tuple (latent variable,output)
-        :param mu: latent variable's predictive mean
-        :param sigma: latent variable's predictive standard deviation
-        ..Note: Only avilable when the output is continuous
-        """
-        assert not self.discrete, "Hessian not available for discrete outputs."
-        cross_derivative = self._d2nlog_mass_dcross(gp=x[0],obs=x[1])
-        return np.array((self._d2nlog_product_dgp2(gp=x[0],obs=x[1],mu=mu,sigma=sigma),cross_derivative,cross_derivative,self._d2nlog_mass_dobs2(obs=x[1],gp=x[0]))).reshape(2,2)
-
-    def _joint_predictive_mode(self,mu,sigma):
-        """
-        Negative logarithm of the joint predictive distribution (latent variable and output).
-
-        :param x: tuple (latent variable,output)
-        :param mu: latent variable's predictive mean
-        :param sigma: latent variable's predictive standard deviation
-        """
-        return sp.optimize.fmin_ncg(self._nlog_joint_predictive_scaled,x0=(mu,self.link.inv_transf(mu)),fprime=self._gradient_nlog_joint_predictive,fhess=self._hessian_nlog_joint_predictive,args=(mu,sigma))
-
-    def predictive_values(self,mu,var,sample=True,sample_size=5000):
-        """
-        Compute  mean, variance and conficence interval (percentiles 5 and 95) of the  prediction
-        :param mu: mean of the latent variable
-        :param var: variance of the latent variable
-        """
-        if isinstance(mu,float) or isinstance(mu,int):
-            mu = [mu]
-            var = [var]
-        pred_mean = []
-        pred_var = []
-        q1 = []
-        q3 = []
-        for m,s in zip(mu,np.sqrt(var)):
-            pred_mean.append(self.predictive_mean(m,s))
-            pred_var.append(self.predictive_variance(m,s,pred_mean[-1]))
-            q1.append(self.predictive_mean(stats.norm.ppf(.025,m,s**2),s))
-            q3.append(self.predictive_mean(stats.norm.ppf(.975,m,s**2),s))
-        pred_mean = np.array(pred_mean)[:,None]
-        pred_var = np.array(pred_var)[:,None]
-        q1 = np.array(q1)[:,None]
-        q3 = np.array(q3)[:,None]
-        return pred_mean, pred_var, q1, q3
diff --git a/GPy/likelihoods/link_functions.py b/GPy/likelihoods/link_functions.py
deleted file mode 100644
index a6434bfb..00000000
--- a/GPy/likelihoods/link_functions.py
+++ /dev/null
@@ -1,100 +0,0 @@
-# Copyright (c) 2012, 2013 Ricardo Andrade
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-from scipy import stats
-import scipy as sp
-import pylab as pb
-from ..util.plot import gpplot
-from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf,inv_std_norm_cdf
-
-class LinkFunction(object):
-    """
-    Link function class for doing non-Gaussian likelihoods approximation
-
-    :param Y: observed output (Nx1 numpy.darray)
-    ..Note:: Y values allowed depend on the likelihood_function used
-    """
-    def __init__(self):
-        pass
-
-class Identity(LinkFunction):
-    """
-    $$
-    g(f) = f
-    $$
-    """
-    def transf(self,mu):
-        return mu
-
-    def inv_transf(self,f):
-        return f
-
-    def dinv_transf_df(self,f):
-        return 1.
-
-    def d2inv_transf_df2(self,f):
-        return 0
-
-
-class Probit(LinkFunction):
-    """
-    $$
-    g(f) = \\Phi^{-1} (mu)
-    $$
-    """
-    def transf(self,mu):
-        return inv_std_norm_cdf(mu)
-
-    def inv_transf(self,f):
-        return std_norm_cdf(f)
-
-    def dinv_transf_df(self,f):
-        return std_norm_pdf(f)
-
-    def d2inv_transf_df2(self,f):
-        return -f * std_norm_pdf(f)
-
-class Log(LinkFunction):
-    """
-    $$
-    g(f) = \log(\mu)
-    $$
-    """
-    def transf(self,mu):
-        return np.log(mu)
-
-    def inv_transf(self,f):
-        return np.exp(f)
-
-    def dinv_transf_df(self,f):
-        return np.exp(f)
-
-    def d2inv_transf_df2(self,f):
-        return np.exp(f)
-
-class Log_ex_1(LinkFunction):
-    """
-    $$
-    g(f) = \log(\exp(\mu) - 1)
-    $$
-    """
-    def transf(self,mu):
-        """
-        function: output space -> latent space
-        """
-        return np.log(np.exp(mu) - 1)
-
-    def inv_transf(self,f):
-        """
-        function: latent space -> output space
-        """
-        return np.log(np.exp(f)+1)
-
-    def dinv_transf_df(self,f):
-        return np.exp(f)/(1.+np.exp(f))
-
-    def d2inv_transf_df2(self,f):
-        aux = np.exp(f)/(1.+np.exp(f))
-        return aux*(1.-aux)
diff --git a/GPy/likelihoods/poisson_likelihood.py b/GPy/likelihoods/poisson_likelihood.py
deleted file mode 100644
index fba89331..00000000
--- a/GPy/likelihoods/poisson_likelihood.py
+++ /dev/null
@@ -1,89 +0,0 @@
-# Copyright (c) 2012, 2013 Ricardo Andrade
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-from scipy import stats,special
-import scipy as sp
-import pylab as pb
-from ..util.plot import gpplot
-from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import link_functions
-from likelihood_functions import LikelihoodFunction
-
-class Poisson(LikelihoodFunction):
-    """
-    Poisson likelihood
-    Y is expected to take values in {0,1,2,...}
-    -----
-    $$
-    L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
-    $$
-    """
-    def __init__(self,link=None):
-        self.discrete = True
-        self.support_limits = (0,np.inf)
-
-        self.analytical_moments = False
-        super(Poisson, self).__init__(link)
-
-    def _mass(self,gp,obs):
-        """
-        Mass (or density) function
-        """
-        return stats.poisson.pmf(obs,self.link.inv_transf(gp))
-
-    def _nlog_mass(self,gp,obs):
-        """
-        Negative logarithm of the un-normalized distribution: factors that are not a function of gp are omitted
-        """
-        return self.link.inv_transf(gp) - obs * np.log(self.link.inv_transf(gp)) + np.log(special.gamma(obs+1))
-
-    #def _preprocess_values(self,Y): #TODO
-
-    def _dnlog_mass_dgp(self,gp,obs):
-        return self.link.dinv_transf_df(gp) * (1. - obs/self.link.inv_transf(gp))
-
-    def _d2nlog_mass_dgp2(self,gp,obs):
-        d2_df = self.link.d2inv_transf_df2(gp)
-        inv_transf = self.link.inv_transf(gp)
-        return obs * ((self.link.dinv_transf_df(gp)/inv_transf)**2 - d2_df/inv_transf) + d2_df
-
-    def _dnlog_mass_dobs(self,obs,gp): #TODO not needed
-        return special.psi(obs+1) -  np.log(self.link.inv_transf(gp))
-
-    def _d2nlog_mass_dobs2(self,obs,gp=None): #TODO not needed
-        return special.polygamma(1,obs)
-
-    def _d2nlog_mass_dcross(self,obs,gp): #TODO not needed
-        return -self.link.dinv_transf_df(gp)/self.link.inv_transf(gp)
-
-    def _mean(self,gp):
-        """
-        Mass (or density) function
-        """
-        return self.link.inv_transf(gp)
-
-    #def _variance(self,gp):
-    #    return self.link.inv_transf(gp)
-
-    def _dmean_dgp(self,gp):
-        return self.link.dinv_transf_df(gp)
-
-    def _d2mean_dgp2(self,gp):
-        return self.link.d2inv_transf_df2(gp)
-
-    def _variance(self,gp):
-        """
-        Mass (or density) function
-        """
-        return self.link.inv_transf(gp)
-
-    #def _variance(self,gp):
-    #    return self.link.inv_transf(gp)
-
-    def _dvariance_dgp(self,gp):
-        return self.link.dinv_transf_df(gp)
-
-    def _d2variance_dgp2(self,gp):
-        return self.link.d2inv_transf_df2(gp)
diff --git a/GPy/models/gp_classification.py b/GPy/models/gp_classification.py
index 2e0d9c4a..d1cf2e00 100644
--- a/GPy/models/gp_classification.py
+++ b/GPy/models/gp_classification.py
@@ -32,7 +32,7 @@ class GPClassification(GP):
 
         if likelihood is None:
             #distribution = GPy.likelihoods.binomial_likelihood.Binomial(link=link)
-            distribution = likelihoods.binomial_likelihood.Binomial()
+            distribution = likelihoods.binomial()
             likelihood = likelihoods.EP(Y, distribution)
         elif Y is not None:
             if not all(Y.flatten() == likelihood.data.flatten()):