docs: def expected_dense_regression_log_prob(A, Sigma, stats):
      """
      Expected log likelihood of p(y | x) where
  
          y ~ N(Ax, Sigma)
  
      and expectation is wrt q(y,x).  We only need expected
      sufficient statistics E[yy.T], E[yx.T], E[xx.T], and n,
      where n is the number of observations.
  
      :param A:      regression matrix
      :param Sigma:  observation covariance
      :param stats:  tuple (E[yy.T], E[yx.T], E[xx.T], n)
      :return:       E[log p(y | x)]
      """
      yyT, yxT, xxT, n = stats[-4:]
  
      contract = 'ij,nij->n' if yyT.ndim == 3 else 'ij,ij->'
      D = A.shape[0]
      Si = np.linalg.inv(Sigma)
      SiA = Si.dot(A)
      ASiA = A.T.dot(SiA)
  
      out = -1. / 2 * np.einsum(contract, ASiA, xxT)
      out += np.einsum(contract, SiA, yxT)
      out += -1. / 2 * np.einsum(contract, Si, yyT)
      out += -D / 2 * np.log(2 * np.pi) + n / 2. * np.linalg.slogdet(Si)[1]
      return out