code: def _expectation(p, mean1, none1, mean2, none2, nghp=None):
      
      return mean1(p.mu)[:, :, None] * mean2(p.mu)[:, None, :]
docs: Compute the expectation:
      expectation[n] = <m1(x_n)^T m2(x_n)>_p(x_n)
          - m1(.), m2(.) :: Constant mean functions
  
      :return: NxQ1xQ2