Fixing bernoulli likelihood for Laplace, fixing Zep for EP, and starting working on quadrature limits

GPy/likelihoods/bernoulli.py

@@ -140,7 +140,7 @@ class Bernoulli(Likelihood):
             Each y_i must be in {0, 1}
         """
         #objective = (inv_link_f**y) * ((1.-inv_link_f)**(1.-y))
-        return np.where(y, inv_link_f, 1.-inv_link_f)
+        return np.where(y==1, inv_link_f, 1.-inv_link_f)
 
     def logpdf_link(self, inv_link_f, y, Y_metadata=None):
         """
@@ -179,7 +179,7 @@ class Bernoulli(Likelihood):
         #grad = (y/inv_link_f) - (1.-y)/(1-inv_link_f)
         #grad = np.where(y, 1./inv_link_f, -1./(1-inv_link_f))
         ff = np.clip(inv_link_f, 1e-9, 1-1e-9)
-        denom = np.where(y, ff, -(1-ff))
+        denom = np.where(y==1, ff, -(1-ff))
         return 1./denom
 
     def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
@@ -205,7 +205,7 @@ class Bernoulli(Likelihood):
         """
         #d2logpdf_dlink2 = -y/(inv_link_f**2) - (1-y)/((1-inv_link_f)**2)
         #d2logpdf_dlink2 = np.where(y, -1./np.square(inv_link_f), -1./np.square(1.-inv_link_f))
-        arg = np.where(y, inv_link_f, 1.-inv_link_f)
+        arg = np.where(y==1, inv_link_f, 1.-inv_link_f)
         ret =  -1./np.square(np.clip(arg, 1e-9, 1e9))
         if np.any(np.isinf(ret)):
             stop
@@ -230,7 +230,7 @@ class Bernoulli(Likelihood):
         #d3logpdf_dlink3 = 2*(y/(inv_link_f**3) - (1-y)/((1-inv_link_f)**3))
         state = np.seterr(divide='ignore')
         # TODO check y \in {0, 1} or {-1, 1}
-        d3logpdf_dlink3 = np.where(y, 2./(inv_link_f**3), -2./((1.-inv_link_f)**3))
+        d3logpdf_dlink3 = np.where(y==1, 2./(inv_link_f**3), -2./((1.-inv_link_f)**3))
         np.seterr(**state)
         return d3logpdf_dlink3
 
@@ -243,8 +243,6 @@ class Bernoulli(Likelihood):
         p = self.predictive_mean(mu, var)
         return [np.asarray(p>(q/100.), dtype=np.int32) for q in quantiles]
 
-
-
     def samples(self, gp, Y_metadata=None):
         """
         Returns a set of samples of observations based on a given value of the latent variable.

GPy/likelihoods/gaussian.py

@@ -67,7 +67,7 @@ class Gaussian(Likelihood):
         """
         return Y
 
-    def _moments_match_ep(self, data_i, tau_i, v_i):
+    def moments_match_ep(self, data_i, tau_i, v_i):
         """
         Moments match of the marginal approximation in EP algorithm
 

GPy/likelihoods/likelihood.py

@@ -49,8 +49,8 @@ class Likelihood(Parameterized):
         """
         return Y.shape[1]
 
-    def _gradients(self,partial):
-        return np.zeros(0)
+    def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None):
+        return np.zeros(self.size)
 
     def update_gradients(self, partial):
         if self.size > 0:
@@ -176,8 +176,10 @@ class Likelihood(Parameterized):
         log_p_ystar = np.array(log_p_ystar).reshape(*y_test.shape)
         return log_p_ystar
 
+    def quad_limits(self):
+        return -np.inf, np.inf
 
-    def _moments_match_ep(self,obs,tau,v):
+    def moments_match_ep(self,obs,tau,v):
         """
         Calculation of moments using quadrature
 
@@ -188,20 +190,27 @@ class Likelihood(Parameterized):
         #Compute first integral for zeroth moment.
         #NOTE constant np.sqrt(2*pi/tau) added at the end of the function
         mu = v/tau
+        sigma2 = 1./tau
+        #Lets do these for now based on the same idea as Gaussian quadrature
+        # i.e. multiply anything by close to zero, and its zero.
+        f_min = mu - 8*np.sqrt(sigma2)
+        f_max = mu + 8*np.sqrt(sigma2)
+
+        # f_min, f_max = self.quad_limits()
         def int_1(f):
             return self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
-        z_scaled, accuracy = quad(int_1, -np.inf, np.inf)
+        z_scaled, accuracy = quad(int_1, f_min, f_max)
 
         #Compute second integral for first moment
         def int_2(f):
             return f*self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
-        mean, accuracy = quad(int_2, -np.inf, np.inf)
+        mean, accuracy = quad(int_2, f_min, f_max)
         mean /= z_scaled
 
         #Compute integral for variance
         def int_3(f):
             return (f**2)*self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
-        Ef2, accuracy = quad(int_3, -np.inf, np.inf)
+        Ef2, accuracy = quad(int_3, f_min, f_max)
         Ef2 /= z_scaled
         variance = Ef2 - mean**2
 

GPy/likelihoods/poisson.py

@@ -28,7 +28,7 @@ class Poisson(Likelihood):
         """
         the expected value of y given a value of f
         """
-        return self.gp_link.transf(gp)
+        return self.gp_link.transf(f)
 
     def pdf_link(self, link_f, y, Y_metadata=None):
         """
@@ -46,7 +46,8 @@ class Poisson(Likelihood):
         :rtype: float
         """
         assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
-        return np.prod(stats.poisson.pmf(y,link_f))
+        return np.exp(self.logpdf_link(link_f, y, Y_metadata))
+        # return np.prod(stats.poisson.pmf(y,link_f))
 
     def logpdf_link(self, link_f, y, Y_metadata=None):
         """