From 361f0a527489f12a3949adc008a650221a455e09 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Thu, 16 Apr 2015 09:25:18 +0100 Subject: [PATCH] Fixed log predictive density, added option for LOO to provide some intemediate variables --- .../latent_function_inference/laplace.py | 18 ++++++++++++------ GPy/likelihoods/likelihood.py | 17 ++++++++++------- 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index ed21f094..aefc82ac 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -40,7 +40,7 @@ class Laplace(LatentFunctionInference): self.first_run = True self._previous_Ki_fhat = None - def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None): + def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None, f_hat=None, W=None, Ki_W_i=None): """ Leave one out log predictive density as found in "Bayesian leave-one-out cross-validation approximations for Gaussian latent variable models" @@ -51,13 +51,19 @@ class Laplace(LatentFunctionInference): if K is None: K = kern.K(X) - f_hat, _ = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) - W = -likelihood.d2logpdf_df2(f_hat, Y, Y_metadata=Y_metadata) + if f_hat is None: + f_hat, _ = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) + + if W is None: + W = -likelihood.d2logpdf_df2(f_hat, Y, Y_metadata=Y_metadata) + + if Ki_W_i is None: + _, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave) + logpdf_dfhat = likelihood.dlogpdf_df(f_hat, Y, Y_metadata=Y_metadata) - K_Wi_i, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave) - - W = np.diagflat(W) + if W.shape[1] == 1: + W = np.diagflat(W) #Eq 14, and 16 var_site = 1./np.diag(W)[:, None] diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 470f5059..34798a35 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -114,21 +114,24 @@ class Likelihood(Parameterized): #Otherwise just pass along None's zipped_values = zip(flat_y_test, flat_mu_star, flat_var_star, [None]*y_test.shape[0]) - def integral_generator(y, m, v, y_m): - """Generate a function which can be integrated to give p(Y*|Y) = int p(Y*|f*)p(f*|Y) df*""" - def f(f_star): + def integral_generator(yi, mi, vi, yi_m): + """Generate a function which can be integrated + to give p(Y*|Y) = int p(Y*|f*)p(f*|Y) df*""" + def f(fi_star): #exponent = np.exp(-(1./(2*v))*np.square(m-f_star)) #from GPy.util.misc import safe_exp #exponent = safe_exp(exponent) #return self.pdf(f_star, y, y_m)*exponent #More stable in the log space - return np.exp(self.logpdf(f_star, y, y_m) -(1./(2*v))*np.square(m-f_star)) + return np.exp(self.logpdf(fi_star, yi, yi_m) + - 0.5*np.log(2*np.pi*vi) + - 0.5*np.square(mi-fi_star)/vi) return f - scaled_p_ystar, accuracy = zip(*[quad(integral_generator(y, m, v, y_m), -np.inf, np.inf) for y, m, v, y_m in zipped_values]) - scaled_p_ystar = np.array(scaled_p_ystar).reshape(-1,1) - p_ystar = scaled_p_ystar/np.sqrt(2*np.pi*var_star) + p_ystar, _ = zip(*[quad(integral_generator(yi, mi, vi, yi_m), -np.inf, np.inf) + for yi, mi, vi, yi_m in zipped_values]) + p_ystar = np.array(p_ystar).reshape(-1, 1) return np.log(p_ystar) def _moments_match_ep(self,obs,tau,v):