Confidence intervals from scratch with Python¶

Libraries and helper functions¶

In [1]:
import math as m

In [2]:
def calc_normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
return (1 + m.erf((x - mu) / m.sqrt(2) / sigma)) / 2

In [3]:
def calc_inverse_normal_cdf(p: float, mu:float = 0, sigma: float = 1, tolerance: float = 1E-5, show_steps=False) -> float:

if p == 0: return -np.inf
if p == 1: return np.inf

# In case it is not a standard normal distribution, calculate the standard normal first and then rescale
if mu != 0 or sigma != 1:
return mu + sigma * calc_inverse_normal_cdf(p, tolerance=tolerance)

low_z = -10
hi_z = 10

if show_steps: print(f"{'':<19}".join(['low_z', 'mid_z', 'hi_z']), "\n")

while hi_z - low_z > tolerance:
mid_z = (low_z + hi_z) / 2
mid_p = calc_normal_cdf(mid_z)

if mid_p < p:
low_z = mid_z
else:
hi_z = mid_z

if show_steps: print("\t".join(map(to_string, [low_z, mid_z, hi_z])))

return mid_z

In [4]:
def normal_upper_bound(probabilty: float, mu: float = 0, sigma: float = 1) -> float:
"""Return z for which P(Z <= z) = probability"""
return calc_inverse_normal_cdf(probabilty, mu, sigma)

def normal_lower_bound(probabilty: float, mu: float = 0, sigma: float = 1) -> float:
"""Return z for which P(Z >= z) = probability"""
return calc_inverse_normal_cdf(1 - probabilty, mu, sigma)

In [5]:
def normal_two_sided_bounds(probability: float, mu: float = 0, sigma: float = 1) -> float:
if probability == 0: return 0, 0

tail_probability = (1 - probability) / 2

lower_bound = normal_upper_bound(tail_probability, mu, sigma)
upper_bound = normal_lower_bound(tail_probability, mu, sigma)

return lower_bound, upper_bound


Example 1¶

We would like to know if a coin is 'fair' or not (i.e. p=0.5)

We flip the coin 1000 times and observe 525 heads.

As we do not know the 'real' probability, we need to rely on our observed results.

In [6]:
p_hat = 525 / 1000
mu = p_hat
sigma = m.sqrt(p_hat * (1 - p_hat) / 1000)

p_hat, mu, sigma

Out[6]:
(0.525, 0.525, 0.015791611697353755)

While we cannot tell if this is really the probabiliy of heads, we can calculate the interval within which we expect the real probability fall with our chosen confidence (here: 95%)

In [7]:
confidence = 0.95
normal_two_sided_bounds(confidence, mu, sigma)

Out[7]:
(0.4940490278129096, 0.5559509721870904)

As 0.5 is within the interval, we cannot reject the hypothesis that the coin is fair.

Example 2¶

In another example we observe 540 heads in 1000 flips.

In [8]:
p_hat = 540 / 1000
mu = p_hat
sigma = m.sqrt(p_hat * (1 - p_hat) / 1000)

p_hat, mu, sigma

Out[8]:
(0.54, 0.54, 0.015760710643876435)

Going through the same steps as above, we get that 0.5 falls outside our 95% confidence interval and therefore we reject the hypothesis that the coin is 'fair'.

In [9]:
normal_two_sided_bounds(confidence, mu, sigma)

Out[9]:
(0.5091095927295919, 0.5708904072704082)