docs: def variational_expectations(self, Fmu, Fvar, Y):
          r"""
          Compute the expected log density of the data, given a Gaussian
          distribution for the function values.
  
          if
              q(f) = N(Fmu, Fvar)
  
          and this object represents
  
              p(y|f)
  
          then this method computes
  
             \int (\log p(y|f)) q(f) df.
  
  
          Here, we implement a default Gauss-Hermite quadrature routine, but some
          likelihoods (Gaussian, Poisson) will implement specific cases.
          """
          return ndiagquad(self.logp,
                           self.num_gauss_hermite_points,
                           Fmu, Fvar, Y=Y)