def _set_gaussian_expected_stats(self, smoothed_mus, smoothed_sigmas, E_xtp1_xtT):
        """
        Both meanfield and VBEM require expected statistics of the continuous latent
        states, x.  This is a helper function to take E[x_t], E[x_t x_t^T] and E[x_{t+1}, x_t^T]
        and compute the expected sufficient statistics for the initial distribution,
        dynamics distribution, and Gaussian observation distribution.
        """
        assert not np.isnan(E_xtp1_xtT).any()
        assert not np.isnan(smoothed_mus).any()
        assert not np.isnan(smoothed_sigmas).any()
        assert smoothed_mus.shape == (self.T, self.D_latent)
        assert smoothed_sigmas.shape == (self.T, self.D_latent, self.D_latent)
        assert E_xtp1_xtT.shape == (self.T-1, self.D_latent, self.D_latent)

        # This is like LDSStates._set_expected_states but doesn't sum over time
        T, D_obs = self.T, self.D_emission
        E_x_xT = smoothed_sigmas + smoothed_mus[:, :, None] * smoothed_mus[:, None, :]
        E_x_uT = smoothed_mus[:, :, None] * self.inputs[:, None, :]
        E_u_uT = self.inputs[:, :, None] * self.inputs[:, None, :]

        E_xu_xuT = np.concatenate((
            np.concatenate((E_x_xT, E_x_uT), axis=2),
            np.concatenate((np.transpose(E_x_uT, (0, 2, 1)), E_u_uT), axis=2)),
            axis=1)
        E_xut_xutT = E_xu_xuT[:-1]
        E_xtp1_xtp1T = E_x_xT[1:]
        E_xtp1_utT = (smoothed_mus[1:, :, None] * self.inputs[:-1, None, :])
        E_xtp1_xutT = np.concatenate((E_xtp1_xtT, E_xtp1_utT), axis=-1)

        # Initial state stats
        self.E_init_stats = (self.smoothed_mus[0], E_x_xT[0], 1.)

        # Dynamics stats
        self.E_dynamics_stats = (E_xtp1_xtp1T, E_xtp1_xutT, E_xut_xutT, np.ones(self.T-1))

        # Emission stats -- special case diagonal noise
        E_yyT = self.data**2 if self.diagonal_noise else self.data[:, :, None] * self.data[:, None, :]
        E_yxT = self.data[:, :, None] * self.smoothed_mus[:, None, :]
        E_yuT = self.data[:, :, None] * self.inputs[:, None, :]
        E_yxuT = np.concatenate((E_yxT, E_yuT), axis=-1)
        E_n = np.ones((T, D_obs)) if self.diagonal_noise else np.ones(T)
        self.E_emission_stats = (E_yyT, E_yxuT, E_xu_xuT, E_n)