docs: def fit(self, obs, pstates, unique_pstates=None):
          """
          Estimate model parameters.
          """
          obs = [np.array(o) for o in obs]
          pstates = [np.array(p) for p in pstates]
  
          # List or array of observation sequences
          assert len(obs) == len(pstates)
          assert obs[0].ndim == 2
          assert pstates[0].ndim == 1
  
          if unique_pstates is not None:
              self._init_pstates(unique_pstates)
          else:
              self._init_pstates(list(set(np.concatenate(pstates))))
  
          # Map the partial states to a unique index
          pstates_idx = [np.array([self.e[p] for p in seq]) for seq in pstates]
  
          if self.init_method == 'rand':
              self._init_random()
          elif self.init_method == 'obs':
              self._init_from_obs(obs, pstates_idx)
  
          self._init_pstate_freqs(pstates_idx)
          self._smooth()
  
          logprob = []
          for i in range(self.max_iter):
              # Expectation step
              stats = self._initialize_sufficient_statistics()
              curr_logprob = 0
              for obs_i, pstates_idx_i in zip(obs, pstates_idx):
                  framelogprob = self._compute_log_likelihood(obs_i, pstates_idx_i)
                  lpr, fwdlattice = self._do_forward_pass(framelogprob, pstates_idx_i)
                  bwdlattice = self._do_backward_pass(framelogprob, pstates_idx_i)
                  gamma = fwdlattice + bwdlattice
                  posteriors = np.exp(gamma.T - logsumexp(gamma, axis=1)).T
                  curr_logprob += lpr
  
                  self._accumulate_sufficient_statistics(stats, obs_i,
                                                         pstates_idx_i, framelogprob, posteriors, fwdlattice, bwdlattice)
  
              logprob.append(curr_logprob)
              self.logprob_ = curr_logprob
  
              # Check for convergence.
              self.n_iter_performed_ = i
              if i > 0:
                  self.logprob_delta = logprob[-1] - logprob[-2]
                  if self.logprob_delta < self.thresh:
                      break
  
              # Maximization step
              self._do_mstep(stats)
  
              # Mix the parameters
              self._smooth()
  
          return self