Determine the joint maximum likelihood estimate

def expected_value(self, feature=None, pstate=None, hstate=None, pstate_prob=None, hstate_prob=None):
        """
        Determine the joint maximum likelihood estimate
        """
        # Use the first feature by default
        if feature is None:
            feature = self.emission_name[0]

        # Will default to marginal pstate if pstate is unknown or None
        pstate_idx = self.e[pstate]

        if pstate is not None and pstate_prob is not None:
            raise Exception('Must provide either pstate or pstate_proba but not both')

        if hstate is not None and hstate_prob is not None:
            raise Exception('Must provide either hstate or hstate_proba but not both')

        # Marginalize pstate using the steady state probas
        if pstate_prob is None:
            pstate_prob = self.pstate_steadyprob

        # Marginalize hstate using the steady state probas
        if hstate_prob is None:
            hstate_prob = self.steadyprob[pstate_idx]

        if pstate is None and hstate is None:
            # Marginalize both pstate and hstate
            w = (pstate_prob[:, np.newaxis] * hstate_prob).flatten()

            if self.emission_name_distr[feature] == 'lognormal':
                return np.sum(w * expected_lognormal(self.emission[feature]['logsigma'].flatten(),
                                                     self.emission[feature]['logmu'].flatten()))
            elif self.emission_name_distr[feature] == 'normal':
                return np.sum(w * expected_normal(self.emission[feature]['mu'].flatten(),
                                                  self.emission[feature]['sigma'].flatten()))
        elif hstate is None:
            # Marginalize hstate
            if self.emission_name_distr[feature] == 'lognormal':
                return np.sum(hstate_prob * expected_lognormal(self.emission[feature]['logsigma'][pstate_idx, :],
                                                               self.emission[feature]['logmu'][pstate_idx, :]))
            elif self.emission_name_distr[feature] == 'normal':
                return np.sum(hstate_prob * expected_normal(self.emission[feature]['mu'][pstate_idx, :],
                                                            self.emission[feature]['sigma'][pstate_idx, :]))

        elif pstate is None:
            # Marginalize pstate
            if self.emission_name_distr[feature] == 'lognormal':
                return np.sum(pstate_prob * expected_lognormal(self.emission[feature]['logsigma'][:, hstate],
                                                               self.emission[feature]['logmu'][:, hstate]))
            elif self.emission_name_distr[feature] == 'normal':
                return np.sum(pstate_prob * expected_normal(self.emission[feature]['mu'][:, hstate],
                                                            self.emission[feature]['sigma'][:, hstate]))
        else:
            if self.emission_name_distr[feature] == 'lognormal':
                return expected_lognormal(self.emission[feature]['logsigma'][pstate_idx, hstate],
                                          self.emission[feature]['logmu'][pstate_idx, hstate])
            elif self.emission_name_distr[feature] == 'normal':
                return expected_normal(self.emission[feature]['mu'][pstate_idx, hstate],
                                       self.emission[feature]['sigma'][pstate_idx, hstate])

        return