Given n tests and k successes, return efficiency and confidence interval.

def confidenceInterval(n, k, alpha = 0.68, errorbar=False):
    """
    Given n tests and k successes, return efficiency and confidence interval.
    """
    try:
        e = float(k) / float(n)
    except ZeroDivisionError:
        return np.nan, [np.nan, np.nan]

    bins = 1000001
    dx = 1. / bins

    efficiency = np.linspace(0, 1, bins)

    # MODIFIED FOR LARGE NUMBERS
    if n + 2 > 1000:
        a = gammalnStirling(n + 2)
    else:
        a = scipy.special.gammaln(n + 2)
    if k + 1 > 1000:
        b = gammalnStirling(k + 1)
    else:
        b = scipy.special.gammaln(k + 1)
    if n - k + 1 > 1000:
        c = gammalnStirling(n - k + 1)
    else:
        c = scipy.special.gammaln(n - k + 1)

    if k == 0:
        p = np.concatenate([[np.exp(a - b - c)],
                               np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])),
                               [0.]])
    elif k == n:
        p = np.concatenate([[0.],
                               np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])),
                               [np.exp(a - b - c)]])
    else:
        p = np.concatenate([[0.],
                               np.exp(a - b - c + (k * np.log(efficiency[1: -1])) + (n - k) * np.log(1. - efficiency[1: -1])),
                               [0.]])

    i = np.argsort(p)[::-1]
    p_i = np.take(p, i)

    s = i[np.cumsum(p_i * dx) < alpha]

    low = min(np.min(s) * dx, e)
    high = max(np.max(s) * dx, e)

    if not errorbar:
        return e, [low, high]
    else:
        return e, [e - low, high - e]