code: def median_interval(data, alpha=_alpha):
      
      q = [100*alpha/2., 50, 100*(1-alpha/2.)]
      lo,med,hi = np.percentile(data,q)
      return interval(med,lo,hi)
docs: Median including bayesian credible interval.