docs: def _quadrature_expectation(p, obj1, feature1, obj2, feature2, num_gauss_hermite_points):
      """
      General handling of quadrature expectations for Gaussians and DiagonalGaussians
      Fallback method for missing analytic expectations
      """
      num_gauss_hermite_points = 100 if num_gauss_hermite_points is None else num_gauss_hermite_points
  
      if obj2 is None:
          eval_func = lambda x: get_eval_func(obj1, feature1)(x)
      elif obj1 is None:
          raise NotImplementedError("First object cannot be None.")
      else:
          eval_func = lambda x: (get_eval_func(obj1, feature1, np.s_[:, :, None])(x) *
                                 get_eval_func(obj2, feature2, np.s_[:, None, :])(x))
  
      if isinstance(p, DiagonalGaussian):
          if isinstance(obj1, kernels.Kernel) and isinstance(obj2, kernels.Kernel) \
                  and obj1.on_separate_dims(obj2):  # no joint expectations required
  
              eKxz1 = quadrature_expectation(p, (obj1, feature1),
                                             num_gauss_hermite_points=num_gauss_hermite_points)
              eKxz2 = quadrature_expectation(p, (obj2, feature2),
                                             num_gauss_hermite_points=num_gauss_hermite_points)
              return eKxz1[:, :, None] * eKxz2[:, None, :]
  
          else:
              cov = tf.matrix_diag(p.cov)
      else:
          cov = p.cov
      return mvnquad(eval_func, p.mu, cov, num_gauss_hermite_points)