docs: 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)