Credible interval of 1 dimensional distributions in Python

credible_interval.py

def credible_interval(xlim, ys, mass=0.95):
    """Given a range, like (0, 20), and the y values for the probability density function,
    returns the set of intervals spanning the minimum range to sum up to the mass."""
    
    indexed = zip(ys, range(len(ys)))
    indexed.sort(reverse=True)
    
    acc = 0
    points = []
    for (y, x) in indexed:
        if acc + y > mass:
            break;
        acc += y
        points.append(x)

    points.sort()

    sections = []
    section = [points[0], points[0]]
    for point in points[1:]:
        if section[1] != (point - 1):
            sections.append(section)
            section = (point, point)
        else:
            section[1] = point
    
    sections.append(section)
    
    xs = np.linspace(xlim[0], xlim[1], len(ys))
    translated = [(xs[lb], xs[ub]) for (lb, ub) in sections]
    return translated