docs: def predict_density(self, Fmu, Fvar, Y, epsilon=None):
          r"""
          Given a Normal distribution for the latent function, and a datum Y,
          compute the log predictive density of Y.
  
          i.e. if
              q(f) = N(Fmu, Fvar)
  
          and this object represents
  
              p(y|f)
  
          then this method computes the predictive density
  
              \log \int p(y=Y|f)q(f) df
  
          Here, we implement a default Monte Carlo routine.
          """
          return self._mc_quadrature(self.logp, Fmu, Fvar, Y=Y, logspace=True, epsilon=epsilon)