def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None, f_hat=None, W=None, Ki_W_i=None):
Ki_f_init = np.zeros_like(Y)<EOL>if K is None:<EOL><INDENT>K = kern.K(X)<EOL><DEDENT>if f_hat is None:<EOL><INDENT>f_hat, _ = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)<EOL><DEDENT>if W is None:<EOL><INDENT>W = -likelihood.d2logpdf_df2(f_hat, Y, Y_metadata=Y_metadata)<EOL><DEDENT>if Ki_W_i is None:<EOL><INDENT>_, _, _, Ki_W_i = self._compute_B_statistics(K, W, likelihood.log_concave)<EOL><DEDENT>logpdf_dfhat = likelihood.dlogpdf_df(f_hat, Y, Y_metadata=Y_metadata)<EOL>if W.shape[<NUM_LIT:1>] == <NUM_LIT:1>:<EOL><INDENT>W = np.diagflat(W)<EOL><DEDENT>var_site = <NUM_LIT:1.>/np.diag(W)[:, None]<EOL>mu_site = f_hat + var_site*logpdf_dfhat<EOL>prec_site = <NUM_LIT:1.>/var_site<EOL>marginal_cov = Ki_W_i<EOL>marginal_mu = marginal_cov.dot(np.diagflat(prec_site)).dot(mu_site)<EOL>marginal_var = np.diag(marginal_cov)[:, None]<EOL>posterior_cav_var = <NUM_LIT:1.>/(<NUM_LIT:1.>/marginal_var - <NUM_LIT:1.>/var_site)<EOL>posterior_cav_mean = posterior_cav_var*((<NUM_LIT:1.>/marginal_var)*marginal_mu - (<NUM_LIT:1.>/var_site)*Y)<EOL>flat_y = Y.flatten()<EOL>flat_mu = posterior_cav_mean.flatten()<EOL>flat_var = posterior_cav_var.flatten()<EOL>if Y_metadata is not None:<EOL><INDENT>Y_metadata_flat = {}<EOL>if Y_metadata is not None:<EOL><INDENT>for key, val in Y_metadata.items():<EOL><INDENT>Y_metadata_flat[key] = np.atleast_1d(val).reshape(-<NUM_LIT:1>, <NUM_LIT:1>)<EOL><DEDENT><DEDENT>zipped_values = []<EOL>for i in range(Y.shape[<NUM_LIT:0>]):<EOL><INDENT>y_m = {}<EOL>for key, val in Y_metadata_flat.items():<EOL><INDENT>if np.isscalar(val) or val.shape[<NUM_LIT:0>] == <NUM_LIT:1>:<EOL><INDENT>y_m[key] = val<EOL><DEDENT>else:<EOL><INDENT>y_m[key] = val[i]<EOL><DEDENT><DEDENT>zipped_values.append((flat_y[i], flat_mu[i], flat_var[i], y_m))<EOL><DEDENT><DEDENT>else:<EOL><INDENT>zipped_values = zip(flat_y, flat_mu, flat_var, [None]*Y.shape[<NUM_LIT:0>])<EOL><DEDENT>def integral_generator(yi, mi, vi, yi_m):<EOL><INDENT>def f(fi_star):<EOL><INDENT>p_fi = np.exp(likelihood.logpdf(fi_star, yi, yi_m)<EOL>- <NUM_LIT:0.5>*np.log(<NUM_LIT:2>*np.pi*vi)<EOL>- <NUM_LIT:0.5>*np.square(mi-fi_star)/vi)<EOL>return p_fi<EOL><DEDENT>return f<EOL><DEDENT>p_ystar, _ = zip(*[quad(integral_generator(y, m, v, yi_m), -np.inf, np.inf)<EOL>for y, m, v, yi_m in zipped_values])<EOL>p_ystar = np.array(p_ystar).reshape(-<NUM_LIT:1>, <NUM_LIT:1>)<EOL>return np.log(p_ystar)<EOL>
Leave one out log predictive density as found in
"Bayesian leave-one-out cross-validation approximations for Gaussian latent variable models"
Vehtari et al. 2014.