code: def predict_density(self, Fmu, Fvar, Y):
          r
          return ndiagquad(self.logp,
                           self.num_gauss_hermite_points,
                           Fmu, Fvar, logspace=True, Y=Y)
docs: 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 Gauss-Hermite quadrature routine, but some
          likelihoods (Gaussian, Poisson) will implement specific cases.