docs: def predict_mean_and_var(self, Fmu, Fvar):
          r"""
          Given a Normal distribution for the latent function,
          return the mean of Y
  
          if
              q(f) = N(Fmu, Fvar)
  
          and this object represents
  
              p(y|f)
  
          then this method computes the predictive mean
  
             \int\int y p(y|f)q(f) df dy
  
          and the predictive variance
  
             \int\int y^2 p(y|f)q(f) df dy  - [ \int\int y p(y|f)q(f) df dy ]^2
  
          Here, we implement a default Gauss-Hermite quadrature routine, but some
          likelihoods (e.g. Gaussian) will implement specific cases.
          """
          integrand2 = lambda *X: self.conditional_variance(*X) + tf.square(self.conditional_mean(*X))
          E_y, E_y2 = ndiagquad([self.conditional_mean, integrand2],
                                self.num_gauss_hermite_points,
                                Fmu, Fvar)
          V_y = E_y2 - tf.square(E_y)
          return E_y, V_y