docs: def find_ML(sampler, modelidx):
      """
      Find Maximum Likelihood parameters as those in the chain with a highest log
      probability.
      """
      index = np.unravel_index(
          np.argmax(sampler.lnprobability), sampler.lnprobability.shape
      )
      MLp = sampler.chain[index]
      if modelidx is not None and hasattr(sampler, "blobs"):
          blob = sampler.blobs[index[1]][index[0]][modelidx]
          if isinstance(blob, u.Quantity):
              modelx = sampler.data["energy"].copy()
              model_ML = blob.copy()
          elif len(blob) == 2:
              modelx = blob[0].copy()
              model_ML = blob[1].copy()
          else:
              raise TypeError("Model {0} has wrong blob format".format(modelidx))
      elif modelidx is not None and hasattr(sampler, "modelfn"):
          blob = _process_blob(
              [sampler.modelfn(MLp, sampler.data)],
              modelidx,
              energy=sampler.data["energy"],
          )
          modelx, model_ML = blob[0], blob[1][0]
      else:
          modelx, model_ML = None, None
  
      MLerr = []
      for dist in sampler.flatchain.T:
          hilo = np.percentile(dist, [16.0, 84.0])
          MLerr.append((hilo[1] - hilo[0]) / 2.0)
      ML = sampler.lnprobability[index]
  
      return ML, MLp, MLerr, (modelx, model_ML)