# Byar's Confidence Interval Calculation

I was looking for advice on confidence intervals. I am currently working on automating the calculation of some of our indicators using Byar's method (outlined in Formula's 4a and 4b on page 7 of this document). I have coded these formulas out in python and I'm getting to-be-expected results with the following functions:

z = 1.96
def byars_lower(count, denominator, rate):
c = 1 / (9 * count)
b = 3 * sqrt(count)
lower_o = count * ((1 - c - (z / b)) ** 3)
lower_ci = (lower_o / denominator) * rate
return lower_ci

def byars_higher(count, denominator, rate):
c = 1 / (9 * (count + 1))
b = 3 * (sqrt(count) + 1)
upper_o = (count + 1) * ((1 - c + (z / b)) ** 3)
upper_ci = (upper_o / denominator) * rate
return upper_ci


The current team practice is to use an excel add-in to generate the CIs that basically have the following formulas:

lower CI = =IF(A2=0,0,IF(A2<389,CHIINV(0.5+95/200,2*A2)/2,A2*(1-1/(9*A2)-NORMSINV(0.5+95/200)/3/SQRT(A2))^3))/F2*C2

upper CI = =IF(A2<389,CHIINV(0.5-95/200,2*A2+2)/2,(A2+1)*(1-1/(9*(A2+1))+NORMSINV(0.5+95/200)/3/SQRT(A2+1))^3)/F2*C2

where: A2 = numerator, F2 = denominator & C2 = rate

There is a small discrepancy between these two methods that I can't reconcile (0.0005% - 0.001%), and they need to match to pass quality assurance. Would it be due to the CHIINV function in the excel formulas that would cause the difference?

So, after some more in depth analysis, there were a few issues:

z = 1.96 is stated in the documentation, but its actual value should have been z = ndtri(1 - alpha / 2) using from scipy.special import ndtri.

Line 3 in the byars_higher function should read b = 3 * (sqrt(count + 1)).

This solves all of the discrepancies for numerators above 389. I had to put a conditional statement in to deal with the numbers under 389. The chi2.ppf function from scipy.stats import chi2 replicates CHIINV from excel. The if statements are:

def byars_lower(count, denominator, rate):
if count < 389:
b = (chi2.ppf((alpha * 2), (count * 2)) / 2)
c = b / denominator
lower_ci = c * rate
return lower_ci

def byars_higher(count, denominator, rate):
if count < 389:
b = chi2.ppf(1 - (alpha / 2), 2 * count + 2) / 2
c = b / denominator
upper_ci = c * rate
return upper_ci


It now replicates the excel function above.