docs: 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