From ab5f3591035c333079499a19a6d35357a9ce7dd0 Mon Sep 17 00:00:00 2001
From: Alan Saul <alan.daniel.saul@gmail.com>
Date: Tue, 14 Apr 2015 16:27:58 +0100
Subject: [PATCH] Changed LOO implementation for Eq 30 instead of 37

---
 .../latent_function_inference/laplace.py      | 22 ++++++++++++++-----
 GPy/likelihoods/likelihood.py                 |  6 ++---
 2 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py
index 19d53505..ed21f094 100644
--- a/GPy/inference/latent_function_inference/laplace.py
+++ b/GPy/inference/latent_function_inference/laplace.py
@@ -57,9 +57,20 @@ class Laplace(LatentFunctionInference):
 
         K_Wi_i, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave)
 
-        #Eq 37
-        posterior_cav_var = 1./(1./np.diag(Ki_W_i) - 1./np.diag(W))[:, None]
-        posterior_cav_mean = f_hat - posterior_cav_var*logpdf_dfhat
+        W = np.diagflat(W)
+
+        #Eq 14, and 16
+        var_site = 1./np.diag(W)[:, None]
+        mu_site = f_hat + var_site*logpdf_dfhat
+        prec_site = 1./var_site
+        #Eq 19
+        marginal_cov = Ki_W_i
+        marginal_mu = marginal_cov.dot(np.diagflat(prec_site)).dot(mu_site)
+        marginal_var = np.diag(marginal_cov)[:, None]
+        #Eq 30 with using site parameters instead of Gaussian site parameters
+        #(var_site instead of sigma^{2} )
+        posterior_cav_var = 1./(1./marginal_var - 1./var_site)
+        posterior_cav_mean = posterior_cav_var*((1./marginal_var)*marginal_mu - (1./var_site)*Y)
 
         flat_y = Y.flatten()
         flat_mu = posterior_cav_mean.flatten()
@@ -90,12 +101,13 @@ class Laplace(LatentFunctionInference):
         def integral_generator(yi, mi, vi, yi_m):
             def f(fi_star):
                 #More stable in the log space
-                return np.exp(likelihood.logpdf(fi_star, yi, yi_m)
+                p_fi = np.exp(likelihood.logpdf(fi_star, yi, yi_m)
                               - 0.5*np.log(2*np.pi*vi)
                               - 0.5*np.square(mi-fi_star)/vi)
+                return p_fi
             return f
 
-        #Eq 25
+        #Eq 30
         p_ystar, _ = zip(*[quad(integral_generator(y, m, v, yi_m), -np.inf, np.inf)
                            for y, m, v, yi_m in zipped_values])
         p_ystar = np.array(p_ystar).reshape(-1, 1)
diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py
index f4b31091..470f5059 100644
--- a/GPy/likelihoods/likelihood.py
+++ b/GPy/likelihoods/likelihood.py
@@ -539,9 +539,9 @@ class Likelihood(Parameterized):
 
         #Parameters are stacked vertically. Must be listed in same order as 'get_param_names'
         # ensure we have gradients for every parameter we want to optimize
-        assert dlogpdf_dtheta.shape[0] == self.size #f, d x num_param array
-        assert dlogpdf_df_dtheta.shape[0] == self.size #f x d x num_param matrix or just f x num_param
-        assert d2logpdf_df2_dtheta.shape[0] == self.size #f x num_param matrix or f x d x num_param matrix, f x f x num_param or f x f x d x num_param
+        assert dlogpdf_dtheta.shape[0] == self.size #num_param array x f, d
+        assert dlogpdf_df_dtheta.shape[0] == self.size #num_param x f x d x matrix or just num_param x f
+        assert d2logpdf_df2_dtheta.shape[0] == self.size #num_param x f matrix or num_param x f x d x matrix, num_param x f x f or num_param x f x f x d
 
         return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta