def compute_distribution_bounds(history, parameter, alpha, run):
    """Returns distribution bounds from pyabc posterior distribution.

    Parameters
    ----------
    history : pyabc.smc
        An object created by :func:`pyabc.abc.run()` or
        :func:`respyabc.respyabc()`.

    parameter : str
        Parameter for which the credible interval should be computed.

    alpha : float
        Level of credibility. Must be between zero and one.

    run : int
        Positive integer determining which pyabc run should be used. If `None`
        last run is used.

    Returns
    -------
    lower: float
        Lower bound of the interval.

    upper: float
        Upper bound of the interval.

    """
    if run is None:
        run = history.max_t

    magnitudes, probabilities = history.get_distribution(m=0, t=run)
    magnitudes["probabilities"] = probabilities
    magnitudes_sorted = magnitudes.sort_values(by=parameter)
    magnitudes_sorted["cum_probabilities"] = magnitudes_sorted["probabilities"].cumsum()
    cut_magnitudes = magnitudes_sorted[
        (magnitudes_sorted["cum_probabilities"] >= alpha / 2)
        & (magnitudes_sorted["cum_probabilities"] <= 1 - alpha / 2)
    ]
    cut_indexed = cut_magnitudes.reset_index(drop=True)
    cut_magnitudes = cut_indexed[parameter]
    lower = cut_magnitudes[0]
    upper = cut_magnitudes[len(cut_magnitudes) - 1]

    return lower, upper