docs: def quadrature_expectation(p, obj1, obj2=None, num_gauss_hermite_points=None):
      """
      Compute the expectation <obj1(x) obj2(x)>_p(x)
      Uses Gauss-Hermite quadrature for approximate integration.
  
      :type p: (mu, cov) tuple or a `ProbabilityDistribution` object
      :type obj1: kernel, mean function, (kernel, features), or None
      :type obj2: kernel, mean function, (kernel, features), or None
      :param int num_gauss_hermite_points: passed to `_quadrature_expectation` to set
                                           the number of Gauss-Hermite points used
      :return: a 1-D, 2-D, or 3-D tensor containing the expectation
      """
      if isinstance(p, tuple):
          assert len(p) == 2
  
          if   p[1].shape.ndims == 2:
              p = DiagonalGaussian(*p)
          elif p[1].shape.ndims == 3:
              p = Gaussian(*p)
          elif p[1].shape.ndims == 4:
              p = MarkovGaussian(*p)
  
      if isinstance(obj1, tuple):
          obj1, feat1 = obj1
      else:
          feat1 = None
  
      if isinstance(obj2, tuple):
          obj2, feat2 = obj2
      else:
          feat2 = None
  
      return _quadrature_expectation(p, obj1, feat1, obj2, feat2, num_gauss_hermite_points)