docs: def _expectation(p, linear_mean, none, kern, feat, nghp=None):
      """
      Compute the expectation:
      expectation[n] = <m(x_n)^T K_{x_n, Z}>_p(x_n)
          - m(x_i) = A x_i + b :: Linear mean function
          - K_{.,.}            :: Kernel function
  
      :return: NxQxM
      """
      with params_as_tensors_for(linear_mean):
          N = p.mu.shape[0].value
          D = p.mu.shape[1].value
          exKxz = expectation(p, mean_functions.Identity(D), (kern, feat), nghp=nghp)
          eKxz = expectation(p, (kern, feat), nghp=nghp)
          eAxKxz = tf.matmul(tf.tile(linear_mean.A[None, :, :], (N, 1, 1)),
                             exKxz, transpose_a=True)
          ebKxz = linear_mean.b[None, :, None] * eKxz[:, None, :]
          return eAxKxz + ebKxz