Quantiles of the binomial and normal distributions

quantiles.py

#!/usr/bin/env python3
# SPDX-License-Identifier: CC0-1.0
from math import exp, factorial, pi, sqrt

def binom(n, k):
	if n - k < 0 or n + k < 0:
		return 0
	assert (n - k) % 2 == 0 and (n + k) % 2 == 0
	return factorial(n) / (2**n * factorial((n-k)//2) * factorial((n+k)//2))
def gauss(n, k):
	# That the normalization actually agrees with de Moivre--Laplace is
	# easiest to verify by looking at the maximum (2 / sqrt(2*pi*n) here,
	# 1 / sqrt(2*pi*n*p*q) there); the exponent is clear from the variance.
	return 2 * exp(-k**2 / (2*n)) / sqrt(2*pi*n)

# <https://en.wikipedia.org/wiki/68-95-99.7_rule>
intervals = [(0.5, 0.382_924_922_548_026),
             (1.0, 0.682_689_492_137_086),
             (1.5, 0.866_385_597_462_284),
             (2.0, 0.954_499_736_103_642),
             (2.5, 0.987_580_669_348_448),
             (3.0, 0.997_300_203_936_740),
             (3.5, 0.999_534_741_841_929),
             (4.0, 0.999_936_657_516_334),
             (4.5, 0.999_993_204_653_751),
             (5.0, 0.999_999_426_696_856)]

n     = 100
ks    = list(range(-n, n+1, 2))
sigma = sqrt(n)

print("SCORE    \tEXACT       \tNORMAL      \tINTEGRAL    ")
for score, integral in intervals:
	exact = sum(binom(n, k) for k in ks if abs(k) <= score * sigma)
	normal = sum(gauss(n, k) for k in ks if abs(k) <= score * sigma)
	print("{:.1f} sigma\t{:.10f}\t{:.10f}\t{:.10f}\t"
	      .format(score, exact, normal, integral))