docs: def state_estimation_ensemble(data, k, n_runs=10, M_list=[], **se_params):
      """
      Runs an ensemble method on the list of M results...
  
      Args:
          data: genes x cells array
          k: number of classes
          n_runs (optional): number of random initializations of state estimation
          M_list (optional): list of M arrays from state estimation
          se_params (optional): optional poisson_estimate_state params
  
      Returns:
          M_new
          W_new
          ll
      """
      if len(M_list)==0:
          M_list = []
          for i in range(n_runs):
              M, W, ll = poisson_estimate_state(data, k, **se_params)
              M_list.append(M)
      M_stacked = np.hstack(M_list)
      M_new, W_new, ll = poisson_estimate_state(M_stacked, k, **se_params)
      W_new = np.dot(data.T, M_new)
      W_new = W_new/W_new.sum(0)
      return M_new, W_new, ll