code: def _hpd_interval(self, x, width):
          
  
          x = np.sort(x)
          n = len(x)
  
          interval_idx_inc = int(np.floor(width * n))
          n_intervals = n - interval_idx_inc
          interval_width = x[interval_idx_inc:] - x[:n_intervals]
  
          if len(interval_width) == 0:
              raise ValueError('Too few elements for interval calculation')
  
          min_idx = np.argmin(interval_width)
          hdi_min = x[min_idx]
          hdi_max = x[min_idx + interval_idx_inc]
  
          index = ['hpd{}_{}'.format(width, x) for x in ['lower', 'upper']]
          return pd.Series([hdi_min, hdi_max], index=index)
docs: Code adapted from pymc3.stats.calc_min_interval:
          https://github.com/pymc-devs/pymc3/blob/master/pymc3/stats.py