# Translate Excel's NORM.S.INV function to Python to calculate confidence interval of an incidence rate

I have an Excel-sheet with some formulas which are able to calculate the confidence interval of an incidence rate. The formulas not from me but they work.

Beside some other stuff the Excel-sheet use the NORM.S.INV function to calculate the incidence rate's confidence interval.

I tried to translate that into Python. My code seems to work because it's output is the same as the Excel-sheet. But I'm not sure what happens here in the details. So I'm not sure if the way is correct or could be improved.

The original Excel-formula (cell names replace by variable names used in the script) is in the comments of that script.

#!/usr/bin/env python3
import math
from statistics import NormalDist

# Number of observed cases
observed_cases_n = 112575
# Size of population
population_n = 752487

# "incidence rate" as cases per 1.000 people
rate_relation = 1000

# Strength of the confidence interval
CI = 95

# Raw incidence rate
ir = observed_cases_n / population_n * rate_relation

# Variance of raw incidence rate
ir_variance = ir * (rate_relation - ir) / population_n

# #########################################################################
# The original Excel formular looks like this.
# I only replaced the cell names with variable names
# I use in that Python script.
# = ir + NORM.S.INV((1 + observed_cases_n / 100) / 2) * SQRT(ir_variance)
# #########################################################################

# Translation of Excel-Function NORM.S.INV()
nd = NormalDist(mu=1, sigma=0.5).inv_cdf((1 + CI / 100) / 2)

# Confidence interval
ci_lower = ir - nd * math.sqrt(ir_variance)
ci_upper = ir + nd * math.sqrt(ir_variance)

print(f'The {CI}% confidence interval is {ci_lower} to {ci_upper}.')


The output is

The 95% confidence interval is 148.78978326403353 to 150.41804358214455.


Are there any suggestions about improvements or mistakes of the mathematical and statistical aspects of that script?

• One good way to test that you have ported a black-box function correctly is to perform a systematic comparison. Since the domain of NORM.S.INV is the open interval $(0,1)$ and this is an infinitely differentiable function, why not simply generate a dataset of values in that interval--perhaps a grid of them supplemented by a few random values--and apply both NORM.S.INV and your solution to them? Its differentiability assures that close agreement on such a dataset will lead to close agreement on all other points in the domain.
– whuber
May 2, 2022 at 14:58