code: def resample_emission_distn(self):
          
          masked_datas = [s.masked_data.tocsc() for s in self.states_list]
          xs = [np.hstack((s.gaussian_states, s.inputs))for s in self.states_list]
  
          for n in range(self.D_obs):
              # Get the nonzero values of the nth column
              rowns = [md.indices[md.indptr[n]:md.indptr[n+1]] for md in masked_datas]
              xns = [x[r] for x,r in zip(xs, rowns)]
              yns = [s.masked_data.getcol(n).data for s in self.states_list]
              maskns = [np.ones_like(y, dtype=bool) for y in yns]
              omegans = [s.omega.getcol(n).data for s in self.states_list]
              self.emission_distn._resample_row_of_emission_matrix(n, xns, yns, maskns, omegans)
docs: Now for the expensive part... the data is stored in a sparse row
          format, which is good for updating the latent states (since we
          primarily rely on dot products with the data, which can be
          efficiently performed for CSR matrices).
  
          However, in order to update the n-th row of the emission matrix,
          we need to know which counts are observed in the n-th column of data.
          This involves converting the data to a sparse column format, which
          can require (time) intensive re-indexing.