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