code: def choose_ancestral_states_mppa(tree, feature, states, force_joint=True):
      lh_feature = get_personalized_feature_name(feature, LH)
      allowed_state_feature = get_personalized_feature_name(feature, ALLOWED_STATES)
      joint_state_feature = get_personalized_feature_name(feature, JOINT_STATE)
      n = len(states)
      _, state2array = get_state2allowed_states(states, False)
      num_scenarios = 1
      unresolved_nodes = 0
      num_states = 0
      for node in tree.traverse():
          marginal_likelihoods = getattr(node, lh_feature)
          marginal_probs = marginal_likelihoods / marginal_likelihoods.sum()
          if force_joint:
              joint_index = getattr(node, joint_state_feature)
              joint_prob = marginal_probs[joint_index]
              marginal_probs = np.hstack((np.sort(np.delete(marginal_probs, joint_index)), [joint_prob]))
          else:
              marginal_probs = np.sort(marginal_probs)
          best_k = n
          best_correstion = np.inf
          for k in range(1, n + 1):
              correction = np.hstack((np.zeros(n - k), np.ones(k) / k)) - marginal_probs
              correction = correction.dot(correction)
              if correction < best_correstion:
                  best_correstion = correction
                  best_k = k
          num_scenarios *= best_k
          num_states += best_k
          if force_joint:
              indices_selected = sorted(range(n),
                                        key=lambda _: (0 if n == joint_index else 1, -marginal_likelihoods[_]))[:best_k]
          else:
              indices_selected = sorted(range(n), key=lambda _: -marginal_likelihoods[_])[:best_k]
          if best_k == 1:
              allowed_states = state2array[indices_selected[0]]
          else:
              allowed_states = np.zeros(len(states), dtype=np.int)
              allowed_states[indices_selected] = 1
              unresolved_nodes += 1
          node.add_feature(allowed_state_feature, allowed_states)
      return num_scenarios, unresolved_nodes, num_states