beginning of adding variational GH quadrature to the likelihood class

GPy/inference/latent_function_inference/laplace.py

@@ -82,7 +82,7 @@ class Laplace(LatentFunctionInference):
 
         #define the objective function (to be maximised)
         def obj(Ki_f, f):
-            return -0.5*np.dot(Ki_f.flatten(), f.flatten()) + likelihood.logpdf(f, Y, Y_metadata=Y_metadata)
+            return -0.5*np.dot(Ki_f.flatten(), f.flatten()) + np.sum(likelihood.logpdf(f, Y, Y_metadata=Y_metadata))
 
         difference = np.inf
         iteration = 0
@@ -152,7 +152,7 @@ class Laplace(LatentFunctionInference):
         Ki_W_i  = K - C.T.dot(C)
 
         #compute the log marginal
-        log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, Y_metadata=Y_metadata) - np.sum(np.log(np.diag(L)))
+        log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, Y_metadata=Y_metadata).sum() - np.sum(np.log(np.diag(L)))
 
         # Compute matrices for derivatives
         dW_df = -likelihood.d3logpdf_df3(f_hat, Y, Y_metadata=Y_metadata) # -d3lik_d3fhat

GPy/likelihoods/bernoulli.py

@@ -136,10 +136,8 @@ class Bernoulli(Likelihood):
         assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape
         #objective = y*np.log(inv_link_f) + (1.-y)*np.log(inv_link_f)
         state = np.seterr(divide='ignore')
-        # TODO check y \in {0, 1} or {-1, 1}
-        objective = np.where(y==1, np.log(inv_link_f), np.log(1-inv_link_f))
-        np.seterr(**state)
-        return np.sum(objective)
+        p = np.where(y==1, inv_link_f, 1.-inv_link_f)
+        return np.log(p)
 
     def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None):
         """

GPy/likelihoods/likelihood.py

@@ -131,6 +131,40 @@ class Likelihood(Parameterized):
 
         return z, mean, variance
 
+    def variational_expectations(self, Y, m, v, gh_points=None):
+        """
+        Use Gauss-Hermite Quadrature to compute 
+
+           E_p(f) [ log p(y|f) ]
+           d/dm E_p(f) [ log p(y|f) ]
+           d/dv E_p(f) [ log p(y|f) ]
+
+        where p(f) is a Gaussian with mean m and variance v. The shapes of Y, m and v should match.
+
+        if no gh_points are passed, we construct them using defualt options
+        """
+
+        if gh_points is None:
+            gh_x, gh_w = np.polynomial.hermite.hermgauss(20)
+
+        shape = m.shape
+        m,v,Y = m.flatten(), v.flatten(), Y.flatten()
+
+        #make a grid of points
+        X = gh_x[None,:]*np.sqrt(2.*v[:,None]) + m[:,None]
+        logp = self.logpdf(X,Y[:,None])
+
+        p = np.clip(p, 1e-9, 1.-1e-9) # for numerical stability
+        N = stats.norm.pdf(X)
+        F = np.log(p).dot(self.gh_w)
+        NoverP = N/p
+        dF_dm = (NoverP*self.Ysign[:,None]).dot(self.gh_w)
+        dF_dv = -0.5*(NoverP**2 + NoverP*X).dot(self.gh_w)
+        return F, dF_dm, dF_dv
+
+
+
+
     def predictive_mean(self, mu, variance, Y_metadata=None):
         """
         Quadrature calculation of the predictive mean: E(Y_star|Y) = E( E(Y_star|f_star, Y) )

GPy/testing/likelihood_tests.py

@@ -353,7 +353,7 @@ class TestNoiseModels(object):
         #print model._get_params()
         np.testing.assert_almost_equal(
                                model.pdf(f.copy(), Y.copy()),
-                               np.exp(model.logpdf(f.copy(), Y.copy()))
+                               np.exp(model.logpdf(f.copy(), Y.copy()).sum())
                                )
 
     @with_setup(setUp, tearDown)