Highest Density Interval

hdi.py

def hdi(frozen_dist, alpha=0.05, **args):
    """ Highest density interval for a unimodel distribution """
    # freeze distribution with given arguments
    def interval_width(hdi_pt_l):
        return frozen_dist.ppf(alpha + hdi_pt_l) - frozen_dist.ppf(hdi_pt_l)

    # find hdi_pt_l that minimizes interval_width
    hdi_pt_l = fmin(interval_width, 1-alpha, ftol=1e-8, disp=False)[0]
    # return interval as array([low, high])
    return frozen_dist.ppf([hdi_pt_l, alpha + hdi_pt_l])