docs: def max_marginal(self, variables=None, evidence=None, elimination_order=None):
          """
          Computes the max-marginal over the variables given the evidence.
  
          Parameters
          ----------
          variables: list
              list of variables over which we want to compute the max-marginal.
          evidence: dict
              a dict key, value pair as {var: state_of_var_observed}
              None if no evidence
          elimination_order: list
              order of variable eliminations (if nothing is provided) order is
              computed automatically
  
          Examples
          --------
          >>> import numpy as np
          >>> import pandas as pd
          >>> from pgmpy.models import BayesianModel
          >>> from pgmpy.inference import VariableElimination
          >>> values = pd.DataFrame(np.random.randint(low=0, high=2, size=(1000, 5)),
          ...                       columns=['A', 'B', 'C', 'D', 'E'])
          >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D'), ('B', 'E')])
          >>> model.fit(values)
          >>> inference = VariableElimination(model)
          >>> phi_query = inference.max_marginal(['A', 'B'])
          """
          if not variables:
              variables = []
          final_distribution = self._variable_elimination(variables, 'maximize',
                                                          evidence=evidence,
                                                          elimination_order=elimination_order)
  
          return np.max(final_distribution.values)