<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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