fraction_uncertainties

fraction_uncertainties.py

def asymmetric_fractional_error_bars(n, total, mass=0.68, centre='peak'):
    """
    returns the error bar for the fraction $x = n/total$ ($f^{+upper}_{-lower}$)
    Use `mass` to specify the errorbar width
    n: number of successes
    total: total number of trials
    :mass: the desired probability mass contained with the error bar
            default=68%, which is ~equivalent to 1 sigma 
    :centre: return either the 'mean', 'median', or most likely of the probability distribution.
             Probably you want 'peak' which is equivalent to n / total 
    :returns: x, (lower, upper) error bars for the fraction x=n/total
    
    example:
    >>> asymmetric_fractional_error_bars(2, 5)  # 2/5 = 40%
    (0.4, (0.15603661866439825, 0.214319177648984))
    >>> asymmetric_fractional_error_bars(np.array([0., 5.]), np.array([5, 5]))  # arrays of 0% and 100%
    (array([0., 1.]),
    (array([0.        , 0.17296289]), array([0.17296289, 0.        ])))
    """
    from scipy.stats import beta
    import numpy as np
    pdf = beta(n + 1, (total - n + 1))
    if centre == 'mean':
        x = pdf.mean()
    elif centre == 'median':
        x = beta.median()
    elif centre == 'peak':
        x = n / total
    else:
        raise ValueError(f"Unknown centre {centre}")
    # when x=0, you still want lower bound to be 0, conserving probability, also when x=1
    low, high = pdf.interval(mass)
    low, high = np.where(x == 0, 0., low), np.where(x == 0, pdf.ppf(mass), high)
    low, high = np.where(x == 1, pdf.ppf(1-mass), low), np.where(x == 1, 1., high)
    return x, (x - low, high - x)