docs: def _expectation(p, kern, none1, none2, none3, nghp=None):
      r"""
      Compute the expectation:
      <\HadamardProd_i diag(Ki_{X[:, active_dims_i], X[:, active_dims_i]})>_p(X)
          - \HadamardProd_i Ki_{.,.} :: Product kernel
          - p                        :: DiagonalGaussian distribution (p.cov NxD)
  
      :return: N
      """
      if not kern.on_separate_dimensions:
          raise NotImplementedError(
              "Product currently needs to be defined on separate dimensions.")  # pragma: no cover
  
      return functools.reduce(tf.multiply, [
          expectation(p, k, nghp=nghp) for k in kern.kernels])