Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Quantiles of the binomial and normal distributions
#!/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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment