Confidence Interval of the Median

median_ci.py

import numpy as np
import scipy.stats as sps
        
    
def median_ci(x, alpha=0.95, method='approximate'):
    """confidence interval of the median"""
    n = len(x)
    y = np.sort(x)
    
    if method[0] == 'a':
        i = sps.norm.ppf(0.5-alpha/2) * np.sqrt(0.25*n) + 0.5 * n
        a = int(np.floor(i))
        b = int(np.ceil(i))
        f = i - a        
        cl = y[a] + (y[b] - y[a]) * f
        cu = y[n-1-b] + (y[n-1-a] - y[n-1-b]) * (1 - f)        
        p_final = alpha
    else:
        p = sps.binom.cdf(np.arange(0, n), n, 0.5)        
        pc = p[::-1] - p    
        i = np.argmin(np.abs(pc-alpha))
        cl, cu = y[i], y[n-i-1]
        p_final = pc[i]
    return cl, cu, p_final