docs: def minimal_imap(self, order):
          """
          Returns a Bayesian Model which is minimal IMap of the Joint Probability Distribution
          considering the order of the variables.
  
          Parameters
          ----------
          order: array-like
              The order of the random variables.
  
          Examples
          --------
          >>> import numpy as np
          >>> from pgmpy.factors.discrete import JointProbabilityDistribution
          >>> prob = JointProbabilityDistribution(['x1', 'x2', 'x3'], [2, 3, 2], np.ones(12)/12)
          >>> bayesian_model = prob.minimal_imap(order=['x2', 'x1', 'x3'])
          >>> bayesian_model
          <pgmpy.models.models.models at 0x7fd7440a9320>
          >>> bayesian_model.edges()
          [('x1', 'x3'), ('x2', 'x3')]
          """
          from pgmpy.models import BayesianModel
  
          def get_subsets(u):
              for r in range(len(u) + 1):
                  for i in itertools.combinations(u, r):
                      yield i
  
          G = BayesianModel()
          for variable_index in range(len(order)):
              u = order[:variable_index]
              for subset in get_subsets(u):
                  if (len(subset) < len(u) and
                          self.check_independence([order[variable_index]], set(u) - set(subset), subset, True)):
                      G.add_edges_from([(variable, order[variable_index]) for variable in subset])
          return G