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