docs: def plot_samples(
      ax,
      sampler,
      modelidx=0,
      sed=True,
      n_samples=100,
      e_unit=u.eV,
      e_range=None,
      e_npoints=100,
      threads=None,
      label=None,
      last_step=False,
  ):
      """Plot a number of samples from the sampler chain.
  
      Parameters
      ----------
      ax : `matplotlib.Axes`
          Axes to plot on.
      sampler : `emcee.EnsembleSampler`
          Sampler
      modelidx : int, optional
          Model index. Default is 0
      sed : bool, optional
          Whether to plot SED or differential spectrum. If `None`, the units of
          the observed spectrum will be used.
      n_samples : int, optional
          Number of samples to plot. Default is 100.
      e_unit : :class:`~astropy.units.Unit` or str parseable to unit
          Unit in which to plot energy axis.
      e_range : list of `~astropy.units.Quantity`, length 2, optional
          Limits in energy for the computation of the model samples and ML model.
          Note that setting this parameter will mean that the samples for the
          model are recomputed and depending on the model speed might be quite
          slow.
      e_npoints : int, optional
          How many points to compute for the model samples and ML model if
          `e_range` is set.
      threads : int, optional
          How many parallel processing threads to use when computing the samples.
          Defaults to the number of available cores.
      last_step : bool, optional
          Whether to only use the positions in the final step of the run (True,
          default) or the whole chain (False).
      """
  
      modelx, model = _read_or_calc_samples(
          sampler,
          modelidx,
          last_step=last_step,
          e_range=e_range,
          e_npoints=e_npoints,
          threads=threads,
      )
      # pick first model sample for units
      f_unit, sedf = sed_conversion(modelx, model[0].unit, sed)
  
      sample_alpha = min(5.0 / n_samples, 0.5)
      for my in model[np.random.randint(len(model), size=n_samples)]:
          ax.loglog(
              modelx.to(e_unit).value,
              (my * sedf).to(f_unit).value,
              color=(0.1,) * 3,
              alpha=sample_alpha,
              lw=1.0,
          )
  
      _plot_MLmodel(ax, sampler, modelidx, e_range, e_npoints, e_unit, sed)
  
      if label is not None:
          ax.set_ylabel(
              "{0} [{1}]".format(label, f_unit.to_string("latex_inline"))
          )