code: def log_predictive_density(self, x_test, y_test, Y_metadata=None):
          
          mu_star, var_star = self._raw_predict(x_test)
          fy = self.warping_function.f(y_test)
          ll_lpd = self.likelihood.log_predictive_density(fy, mu_star, var_star, Y_metadata=Y_metadata)
          return ll_lpd + np.log(self.warping_function.fgrad_y(y_test))
docs: Calculation of the log predictive density. Notice we add
          the jacobian of the warping function here.
  
          .. math:
              p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})
  
          :param x_test: test locations (x_{*})
          :type x_test: (Nx1) array
          :param y_test: test observations (y_{*})
          :type y_test: (Nx1) array
          :param Y_metadata: metadata associated with the test points