massive changes

GPy/likelihoods/ep.py

@@ -24,9 +24,18 @@ class EP(likelihood):
 
         #Initial values - Likelihood approximation parameters:
         #p(y|f) = t(f|tau_tilde,v_tilde)
+        #TODO restore
         self.tau_tilde = np.zeros(self.N)
         self.v_tilde = np.zeros(self.N)
 
+        #_gp = self.LikelihoodFunction.link.transf(self.data)
+        #_mean = self.LikelihoodFunction._mean(_gp)
+        #_variance = self.LikelihoodFunction._variance(_gp)
+        #self.tau_tilde = 1./_variance
+        #self.tau_tilde[_variance== 0] = 1.
+        #self.v_tilde = _mean*self.tau_tilde
+
+
         #initial values for the GP variables
         self.Y = np.zeros((self.N,1))
         self.covariance_matrix = np.eye(self.N)
@@ -38,16 +47,17 @@ class EP(likelihood):
         self.trYYT = 0.
 
     def restart(self):
+        #FIXME
         self.tau_tilde = np.zeros(self.N)
         self.v_tilde = np.zeros(self.N)
-        self.Y = np.zeros((self.N,1))
+        #self.Y = np.zeros((self.N,1))
-        self.covariance_matrix = np.eye(self.N)
+        #self.covariance_matrix = np.eye(self.N)
-        self.precision = np.ones(self.N)[:,None]
+        #self.precision = np.ones(self.N)[:,None]
-        self.Z = 0
+        #self.Z = 0
-        self.YYT = None
+        #self.YYT = None
-        self.V = self.precision * self.Y
+        #self.V = self.precision * self.Y
-        self.VVT_factor = self.V
+        #self.VVT_factor = self.V
-        self.trYYT = 0.
+        #self.trYYT = 0.
 
     def predictive_values(self,mu,var,full_cov):
         if full_cov:
@@ -78,6 +88,8 @@ class EP(likelihood):
         self.VVT_factor = self.V
         self.trYYT = np.trace(self.YYT)
 
+        #a = kjkjkjkj
+
     def fit_full(self,K):
         """
         The expectation-propagation algorithm.
@@ -117,15 +129,103 @@ class EP(likelihood):
                 self.v_[i] = mu[i]/Sigma[i,i] - self.eta*self.v_tilde[i]
                 #Marginal moments
                 self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.LikelihoodFunction.moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
+
+
+                #DELETE
+                """
+                import pylab as pb
+                from scipy import stats
+                import scipy as sp
+                import link_functions
+                from constructors import *
+
+                link = link_functions.Log_ex_1()
+                distribution = poisson(link=link)
+                gp = np.linspace(-3,50,100)
+                #distribution = binomial()
+                #gp = np.linspace(-3,3,100)
+
+                y = self._transf_data[i]
+                tau_ = self.tau_[i]
+                v_ = self.v_[i]
+                sigma2_ = np.sqrt(1./tau_)
+                mu_ = v_/tau_
+
+                gaussian = stats.norm.pdf(gp,loc=mu_,scale=np.sqrt(sigma2_))
+                non_gaussian = np.array([distribution._mass(gp_i,y) for gp_i in gp])
+                prod = np.array([distribution._product(gp_i,y,mu_,np.sqrt(sigma2_)) for gp_i in gp])
+                my_Z_hat,my_mu_hat,my_sigma2_hat = distribution.moments_match(y,tau_,v_)
+                proxy = stats.norm.pdf(gp,loc=my_mu_hat,scale=np.sqrt(my_sigma2_hat))
+
+
+                new_sigma2_tilde = 1./self.tau_tilde[i]
+                new_mu_tilde = self.v_tilde[i]/self.tau_tilde[i]
+                new_Z_tilde = self.Z_hat[i]*np.sqrt(2*np.pi)*np.sqrt(sigma2_+new_sigma2_tilde)*np.exp(.5*(mu_-new_mu_tilde)**2/(sigma2_+new_sigma2_tilde))
+                bad_gaussian = stats.norm.pdf(gp,self.v_tilde[i]/self.tau_tilde[i],np.sqrt(1./self.tau_tilde[i]))
+                new_gaussian = stats.norm.pdf(gp,new_mu_tilde,np.sqrt(new_sigma2_tilde))*new_Z_tilde
+                #new_gaussian = stats.norm.pdf(gp,_mu_tilde,np.sqrt(_sigma2_tilde))*_Z_tilde
+
+                _sigma2_tilde = 1./(1./(my_sigma2_hat) - 1./sigma2_)
+                _mu_tilde = (my_mu_hat/my_sigma2_hat - mu_/sigma2_)*_sigma2_tilde
+                _Z_tilde = my_Z_hat*np.sqrt(2*np.pi)*np.sqrt(sigma2_+_sigma2_tilde)*np.exp(.5*(mu_ - _mu_tilde)**2/(sigma2_ + _sigma2_tilde))
+
+                fig1 = pb.figure(figsize=(15,5))
+                ax1 = fig1.add_subplot(131)
+                ax1.grid(True)
+                #pb.plot(gp,bad_gaussian,'b--',linewidth=1.5)
+                #pb.plot(gp,non_gaussian,'b-',linewidth=1.5)
+                pb.plot(gp,new_gaussian,'r--',linewidth=1.5)
+                pb.title('Likelihood: $p(y_i|f_i)$',fontsize=22)
+
+                ax2 = fig1.add_subplot(132)
+                ax2.grid(True)
+                pb.plot(gp,gaussian,'b-',linewidth=1.5)
+                pb.title('Cavity distribution: $q_{-i}(f_i)$',fontsize=22)
+
+                ax3 = fig1.add_subplot(133)
+                ax3.grid(True)
+                pb.plot(gp,prod,'b--',linewidth=1.5)
+
+                pb.plot(gp,proxy*my_Z_hat,'r-',linewidth=1.5)
+
+                pb.title('Approximation: $\mathcal{N}(f_i|\hat{\mu}_i,\hat{\sigma}_i^2) \hat{Z}_i$',fontsize=22)
+                pb.legend(('Exact','Approximation'),frameon=False)
+
+                print 'i',i
+                print 'v/tau _tilde', self.v_tilde[i], self.tau_tilde[i]
+                print 'v/tau _', self.v_[i], self.tau_[i]
+                print 'Z/mu/sigma2 _hat', self.Z_hat[i], mu_hat[i], sigma2_hat[i]
+                pb.plot(gp,new_gaussian*gaussian,'k-')
+
+                a = kj
+                break
+                """
+                #DELETE
+
+
+
+
                 #Site parameters update
-                Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i])
+                Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i]) #FIXME
-                Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i])
+                Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i]) #FIXME
                 self.tau_tilde[i] += Delta_tau
                 self.v_tilde[i] += Delta_v
+
+                #new_tau = self.delta/self.eta*(1./sigma2_hat[i] - self.tau_[i])
+                #new_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - self.v_[i])
+                #Delta_tau = new_tau - self.tau_tilde[i]
+                #Delta_v = new_v - self.v_tilde[i]
+                #self.tau_tilde[i] += Delta_tau
+                #self.v_tilde[i] += Delta_v
+
                 #Posterior distribution parameters update
                 DSYR(Sigma,Sigma[:,i].copy(), -float(Delta_tau/(1.+ Delta_tau*Sigma[i,i])))
                 mu = np.dot(Sigma,self.v_tilde)
                 self.iterations += 1
+
+
+
+
             #Sigma recomptutation with Cholesky decompositon
             Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K
             B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K
@@ -138,6 +238,11 @@ class EP(likelihood):
             self.np1.append(self.tau_tilde.copy())
             self.np2.append(self.v_tilde.copy())
 
+            ##DELETE
+            #pb.vlines(mu[i],0,max(prod))
+            #break
+            #DELETE
+
         return self._compute_GP_variables()
 
     def fit_DTC(self, Kmm, Kmn):

GPy/likelihoods/noise_models/__init__.py

@@ -0,0 +1,4 @@
+import likelihood_functions
+import binomial_likelihood
+import poisson_likelihood
+import link_functions

GPy/likelihoods/noise_models/binomial_likelihood.py

@@ -0,0 +1,111 @@
+# 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
+from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
+import link_functions
+from likelihood_functions import NoiseModel
+
+class Binomial(NoiseModel):
+    """
+    Probit likelihood
+    Y is expected to take values in {-1,1}
+    -----
+    $$
+    L(x) = \\Phi (Y_i*f_i)
+    $$
+    """
+    def __init__(self,link=None,analytical_moments=False):
+        super(Binomial, self).__init__(link,analytical_moments)
+
+    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.link.d2inv_transf_df2(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
+    """

GPy/likelihoods/noise_models/likelihood_functions.py

@@ -0,0 +1,349 @@
+# 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 GPy.util.plot import gpplot
+from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
+import link_functions
+
+
+class NoiseModel(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,analytical_moments=False):
+        #assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME
+        self.link = link
+        self.analytical_moments = analytical_moments
+        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.vstack(pred_mean)
+        pred_var = np.vstack(pred_var)
+        q1 = np.vstack(q1)
+        q3 = np.vstack(q3)
+        return pred_mean, pred_var, q1, q3

GPy/likelihoods/noise_models/link_functions.py

@@ -0,0 +1,99 @@
+# 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 GPy.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(1.+np.exp(f))
+
+    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)

GPy/likelihoods/noise_models/poisson_likelihood.py

@@ -0,0 +1,92 @@
+# 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 GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
+import link_functions
+from likelihood_functions import NoiseModel
+
+class Poisson(NoiseModel):
+    """
+    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,analytical_moments=False):
+        #self.discrete = True
+        #self.support_limits = (0,np.inf)
+
+        #self.analytical_moments = False
+        super(Poisson, self).__init__(link,analytical_moments)
+
+    def _preprocess_values(self,Y): #TODO
+        self.scale = .5*Y.max()
+        self.shift = Y.mean()
+        return (Y - self.shift)/self.scale
+
+    def _mass(self,gp,obs):
+        """
+        Mass (or density) function
+        """
+        obs = obs*self.scale + self.shift
+        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 _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)