hdi.py

def _highest_density_interval(self, pmf, p=.9, title=''):
    # If we pass a DataFrame, just call this recursively on the columns
    if(isinstance(pmf, pd.DataFrame)):
        return pd.DataFrame([self._highest_density_interval(pmf[col], title=str(col)) for col in pmf],
                            index=pmf.columns)

    # Broadcast the probability distribution to an artificial set of samples by
    # repeating each index value N times where N = probability sample_precision
    sample_precision = 1000000
    samples_repeats = np.array(pmf.values * sample_precision)#.astype(int)
    samples = np.repeat(pmf.index, np.array(pmf.values * sample_precision).astype(int))

    # Get HDI
    hdi = pymc3.stats.hpd(samples, credible_interval=p)

    low = hdi[0]
    high = hdi[1]

    return pd.Series([low, high], index=['Rt_low', 'Rt_high'])