Calculation of the log predictive density

        .. math:
            p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})

        :param y_test: test observations (y_{*})
        :type y_test: (Nx1) array
        :param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
        :type mu_star: (Nx1) array
        :param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
        :type var_star: (Nx1) array

def log_predictive_density(self, y_test, mu_star, var_star, Y_metadata=None):
        """
        Calculation of the log predictive density

        .. math:
            p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})

        :param y_test: test observations (y_{*})
        :type y_test: (Nx1) array
        :param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
        :type mu_star: (Nx1) array
        :param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
        :type var_star: (Nx1) array
        """
        assert y_test.shape==mu_star.shape
        assert y_test.shape==var_star.shape
        assert y_test.shape[1] == 1

        flat_y_test = y_test.flatten()
        flat_mu_star = mu_star.flatten()
        flat_var_star = var_star.flatten()

        if Y_metadata is not None:
            #Need to zip individual elements of Y_metadata aswell
            Y_metadata_flat = {}
            if Y_metadata is not None:
                for key, val in Y_metadata.items():
                    Y_metadata_flat[key] = np.atleast_1d(val).reshape(-1,1)

            zipped_values = []

            for i in range(y_test.shape[0]):
                y_m = {}
                for key, val in Y_metadata_flat.items():
                    if np.isscalar(val) or val.shape[0] == 1:
                        y_m[key] = val
                    else:
                        #Won't broadcast yet
                        y_m[key] = val[i]
                zipped_values.append((flat_y_test[i], flat_mu_star[i], flat_var_star[i], y_m))
        else:
            #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(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*vi))*np.square(mi-fi_star))
                #from GPy.util.misc import safe_exp
                #exponent = safe_exp(exponent)
                #res = safe_exp(self.logpdf(fi_star, yi, yi_m))*exponent

                #More stable in the log space
                res = np.exp(self.logpdf(fi_star, yi, yi_m)
                              - 0.5*np.log(2*np.pi*vi)
                              - 0.5*np.square(fi_star-mi)/vi)
                if not np.isfinite(res):
                    import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
                return res

            return f

        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(*y_test.shape)
        return np.log(p_ystar)