docs: def expected_hmm_logprob(pi_0, trans_matrix, stats):
      """
      :param pi_0:          initial distribution
      :param trans_matrix:  transition matrix
      :param stats:         tuple (E[z_t], \sum_t E[z_t z_{t+1}.T])
  
      :return:  E_{q(z)} [ log p(z) ]
      """
      E_z, sum_E_ztztp1T, _ = stats
      T, K = E_z.shape
      assert sum_E_ztztp1T.shape == (K, K)
  
      out = 0
      out += np.dot(E_z[0], np.log(pi_0))
      out += np.sum(sum_E_ztztp1T * np.log(trans_matrix))
      return out