docs: def plot_distribution(samples, label, figure=None):
      """ Plot a distribution and print statistics about it"""
  
      from scipy import stats
      import matplotlib.pyplot as plt
  
      quant = [16, 50, 84]
      quantiles = dict(six.moves.zip(quant, np.percentile(samples, quant)))
      std = np.std(samples)
  
      if isinstance(samples[0], u.Quantity):
          unit = samples[0].unit
      else:
          unit = ""
  
      if isinstance(std, u.Quantity):
          std = std.value
  
      dist_props = "{label} distribution properties:\n \
      $-$ median: ${median}$ {unit}, std: ${std}$ {unit}\n \
      $-$ Median with uncertainties based on \n \
        the 16th and 84th percentiles ($\sim$1$\sigma$):\n\
            {label} = {value_error} {unit}".format(
          label=label,
          median=_latex_float(quantiles[50]),
          std=_latex_float(std),
          value_error=_latex_value_error(
              quantiles[50],
              quantiles[50] - quantiles[16],
              quantiles[84] - quantiles[50],
          ),
          unit=unit,
      )
  
      if figure is None:
          f = plt.figure()
      else:
          f = figure
  
      ax = f.add_subplot(111)
      f.subplots_adjust(bottom=0.40, top=0.93, left=0.06, right=0.95)
  
      f.text(0.2, 0.27, dist_props, ha="left", va="top")
  
      histnbins = min(max(25, int(len(samples) / 100.0)), 100)
      xlabel = "" if label is None else label
  
      if isinstance(samples, u.Quantity):
          samples_nounit = samples.value
      else:
          samples_nounit = samples
  
      n, x, _ = ax.hist(
          samples_nounit,
          histnbins,
          histtype="stepfilled",
          color=color_cycle[0],
          lw=0,
          density=True,
      )
  
      kde = stats.kde.gaussian_kde(samples_nounit)
      ax.plot(x, kde(x), color="k", label="KDE")
  
      ax.axvline(
          quantiles[50],
          ls="--",
          color="k",
          alpha=0.5,
          lw=2,
          label="50% quantile",
      )
      ax.axvspan(
          quantiles[16],
          quantiles[84],
          color=(0.5,) * 3,
          alpha=0.25,
          label="68% CI",
          lw=0,
      )
      # ax.legend()
      for l in ax.get_xticklabels():
          l.set_rotation(45)
      # [l.set_rotation(45) for l in ax.get_yticklabels()]
      if unit != "":
          xlabel += " [{0}]".format(unit)
      ax.set_xlabel(xlabel)
      ax.set_title("posterior distribution of {0}".format(label))
      ax.set_ylim(top=n.max() * 1.05)
  
      return f