code: def vbem_aBl(self):
          
          vbem_aBl = np.zeros((self.T, self.num_states))
          ids, dds, eds = self.init_dynamics_distns, self.dynamics_distns, self.emission_distns
  
          for k, (id, dd) in enumerate(zip(ids, dds)):
              vbem_aBl[0, k] = expected_gaussian_logprob(id.mu, id.sigma, self.E_init_stats)
              vbem_aBl[:-1, k] += expected_regression_log_prob(dd, self.E_dynamics_stats)
  
          if self.single_emission:
              ed = self.emission_distns[0]
              vbem_aBl += expected_regression_log_prob(ed, self.E_emission_stats)[:,None]
          else:
              for k, ed in enumerate(self.emission_distns):
                  vbem_aBl[:, k] += expected_regression_log_prob(ed, self.E_emission_stats)
  
          vbem_aBl[np.isnan(vbem_aBl).any(1)] = 0.
          return vbem_aBl
docs: These are the expected log likelihoods (node potentials)
          as seen from the discrete states.  In other words,
          E_{q(x)} [log p(y, x | z)]