docs: def _quadrature_expectation(p, obj1, feature1, obj2, feature2, num_gauss_hermite_points):
      """
      Handling of quadrature expectations for Markov Gaussians (useful for time series)
      Fallback method for missing analytic expectations wrt Markov Gaussians
      Nota Bene: obj1 is always associated with x_n, whereas obj2 always with x_{n+1}
                 if one requires e.g. <x_{n+1} K_{x_n, Z}>_p(x_{n:n+1}), compute the
                 transpose and then transpose the result of the expectation
      """
      num_gauss_hermite_points = 40 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)
          mu, cov = p.mu[:-1], p.cov[0, :-1]  # cross covariances are not needed
      elif obj1 is None:
          eval_func = lambda x: get_eval_func(obj2, feature2)(x)
          mu, cov = p.mu[1:], p.cov[0, 1:]  # cross covariances are not needed
      else:
          eval_func = lambda x: (get_eval_func(obj1, feature1, np.s_[:, :, None])(tf.split(x, 2, 1)[0]) *
                                 get_eval_func(obj2, feature2, np.s_[:, None, :])(tf.split(x, 2, 1)[1]))
          mu = tf.concat((p.mu[:-1, :], p.mu[1:, :]), 1)  # Nx2D
          cov_top = tf.concat((p.cov[0, :-1, :, :], p.cov[1, :-1, :, :]), 2)  # NxDx2D
          cov_bottom = tf.concat((tf.matrix_transpose(p.cov[1, :-1, :, :]), p.cov[0, 1:, :, :]), 2)
          cov = tf.concat((cov_top, cov_bottom), 1)  # Nx2Dx2D
  
      return mvnquad(eval_func, mu, cov, num_gauss_hermite_points)