docs: def query(self, variables, evidence=None, elimination_order=None, joint=True):
          """
          Parameters
          ----------
          variables: list
              list of variables for which you want to compute the probability
  
          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
  
          joint: boolean (default: True)
              If True, returns a Joint Distribution over `variables`.
              If False, returns a dict of distributions over each of the `variables`.
  
          Examples
          --------
          >>> from pgmpy.inference import VariableElimination
          >>> from pgmpy.models import BayesianModel
          >>> import numpy as np
          >>> import pandas as pd
          >>> 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.query(['A', 'B'])
          """
          return self._variable_elimination(variables, 'marginalize',
                                            evidence=evidence, elimination_order=elimination_order,
                                            joint=joint)