further corrections

GPy/likelihoods/binomial_likelihood.py

@@ -23,17 +23,15 @@ class Binomial(LikelihoodFunction):
     def __init__(self,link=None):
         self.discrete = True
         self.support_limits = (0,1)
-        self._analytical = link_functions.Probit
+
         if not link:
-            link = self._analytical
+            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 _mass(self,gp,obs):
-        pass
-
-    def _nlog_mass(self,gp,obs):
-        pass
-
     def _preprocess_values(self,Y):
         """
         Check if the values of the observations correspond to the values
@@ -66,12 +64,51 @@ class Binomial(LikelihoodFunction):
     def _predictive_mean_analytical(self,mu,sigma):
         return stats.norm.cdf(mu/np.sqrt(1+sigma**2))
 
-    def predictive_values(self,mu,var):
+    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):
         """
-        Compute  mean, variance and conficence interval (percentiles 5 and 95) of the  prediction
+        Mass (or density) function
-        :param mu: mean of the latent variable
-        :param var: variance of the latent variable
         """
+        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))
@@ -83,3 +120,4 @@ class Binomial(LikelihoodFunction):
         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/likelihood_functions.py

@@ -19,12 +19,13 @@ class LikelihoodFunction(object):
     ..Note:: Y values allowed depend on the LikelihoodFunction used
     """
     def __init__(self,link):
-        if link == self._analytical:
+        #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:
-            assert isinstance(link,link_functions.LinkFunction)
-            self.link = link
             self.moments_match = self._moments_match_numerical
             self.predictive_mean = self._predictive_mean_numerical
 
@@ -258,6 +259,7 @@ class LikelihoodFunction(object):
         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])
@@ -267,6 +269,7 @@ class LikelihoodFunction(object):
         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)
@@ -323,6 +326,8 @@ class LikelihoodFunction(object):
     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]

GPy/likelihoods/poisson_likelihood.py

@@ -23,9 +23,8 @@ class Poisson(LikelihoodFunction):
     def __init__(self,link=None):
         self.discrete = True
         self.support_limits = (0,np.inf)
-        self._analytical = None
+
-        if not link:
+        self.analytical_moments = False
-            link = link_functions.Log()
         super(Poisson, self).__init__(link)
 
     def _mass(self,gp,obs):
@@ -34,15 +33,14 @@ class Poisson(LikelihoodFunction):
         """
         return stats.poisson.pmf(obs,self.link.inv_transf(gp))
 
-    def _percentile(self,x,gp,*args): #TODO *args
-        return stats.poisson.ppf(x,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))
 
@@ -66,8 +64,8 @@ class Poisson(LikelihoodFunction):
         """
         return self.link.inv_transf(gp)
 
-    def _variance(self,gp):
+    #def _variance(self,gp):
-        return self.link.inv_transf(gp)
+    #    return self.link.inv_transf(gp)
 
     def _dmean_dgp(self,gp):
         return self.link.dinv_transf_df(gp)
@@ -81,8 +79,8 @@ class Poisson(LikelihoodFunction):
         """
         return self.link.inv_transf(gp)
 
-    def _variance(self,gp):
+    #def _variance(self,gp):
-        return self.link.inv_transf(gp)
+    #    return self.link.inv_transf(gp)
 
     def _dvariance_dgp(self,gp):
         return self.link.dinv_transf_df(gp)

GPy/models/gp_classification.py

@@ -31,7 +31,8 @@ class GPClassification(GP):
             kernel = kern.rbf(X.shape[1])
 
         if likelihood is None:
-            distribution = likelihoods.likelihood_functions.Binomial()
+            #distribution = GPy.likelihoods.binomial_likelihood.Binomial(link=link)
+            distribution = likelihoods.binomial_likelihood.Binomial()
             likelihood = likelihoods.EP(Y, distribution)
         elif Y is not None:
             if not all(Y.flatten() == likelihood.data.flatten()):