# SPDX-License-Identifier: MIT # statistics — decide what data shows, and calculate statistics and # probability. # # Two halves that share nothing but the subject. # # STRUCTURE. You supply three kinds of record as facts. # # Strata of a study (`statistics.Stratum`): for each stratum, the successes and # the trials of two arms, A and B. The study (`statistics.Study`) is folded # from its strata. The rules conclude each arm's rate in each stratum and in # the pooled study, and `paradox` where the pooled comparison reverses what # every stratum shows (Simpson's paradox): # # import statistics # statistics.stratum("stones", "small", 81, 87, 234, 270) # statistics.stratum("stones", "large", 192, 263, 55, 80) # print(statistics.paradoxes(), statistics.pooled("stones")) # # A screening test (`Screen`) by the prevalence of the condition and the # test's sensitivity and specificity. The rules conclude the predictive value # of a positive result and of a negative one, and whether a positive result is # more likely wrong than right: # # statistics.screen("s1", 0.01, 0.9, 0.9) # print(statistics.predictive_values("s1")) # # A significance test (`Significance`) by its p-value and the level it is held # to. The rules conclude whether the result is significant at that level: # # statistics.significance("t1", 0.049, 0.05) # print(statistics.significant_tests()) # # Each conclusion answers `why`, and `why_not` names the premise that failed. # A corrected count withdraws what rested on it: a stratum rewritten can # remove a paradox. Where two sources disagree about a fact, both values are # reasoned about, and `unsettled()` names what they disagree about. # # The pooled counts are the exception, since a sum adds every value it reads: # a stratum whose successes a second source writes, even the same number, # would be counted twice, 354 successes in 350 trials, and a paradox both # sources show would vanish. So while a stratum is bound to its study twice by # one of its four counts, as a second value of that count binds it, that count # is disputed and not pooled, and neither is what rests on it, the rate and # the paradox; a count both sources agree on is still pooled. # `disputed_studies()` names the study, `pooled()` refuses it, and `why_not` on # the pooled count names the dispute. Forgetting the claim that is not to # stand settles it. A stratum whose study two sources disagree on counts in each # study a value names, so each of those pools is one source's. # # CALCULATION. Value scripts that return a value and write nothing: the mean, # the sample variance by Welford's method and the median; least squares and # Pearson's r; the normal distribution through `math.erfc` and its inverse by # bisection; the binomial and Poisson probabilities; the birthday problem; # Welch's t; the chi-square statistic and its p for one and two degrees of # freedom; Bayes' rule; and a seeded Lehmer generator with a Monte Carlo # estimate of pi. DKE Python has no `random` module, so the generator is the # module's own, and a run with the same seed repeats exactly. # # A sample variance needs two values and a rate needs a trial: a stratum or a # list that is too small is the caller's to avoid, and divides by zero. # # Sources, fetched 2026-11-08: # - The kidney-stone counts in the example above are those of Charig et al., # Br Med J 1986; 292(6524): 879-882, as given at # https://en.wikipedia.org/wiki/Simpson%27s_paradox, which also says that # Simpson described the paradox in 1951. # - Anscombe's quartet is Francis Anscombe's of 1973, "Graphs in Statistical # Analysis", The American Statistician 27(1): 17-21; all four sets have # mean x 9, mean y 7.50, line y = 3.00 + 0.500x and r 0.816 # (https://en.wikipedia.org/wiki/Anscombe%27s_quartet). # - Gosset developed the t-test at the Guinness brewery and published it in # Biometrika in 1908 as "Student" # (https://en.wikipedia.org/wiki/Student%27s_t-test). # - 23 people give a chance of about 50.7% that two share a birthday # (https://en.wikipedia.org/wiki/Birthday_problem). # - The Lehmer generator's modulus is 2^31 - 1; Park and Miller (1988) gave # the multiplier 16807 and Park, Miller and Stockmeyer (1993) suggested # 48271, which this module uses # (https://en.wikipedia.org/wiki/Lehmer_random_number_generator). # # Floats are IEEE doubles. `math.erf`, `math.erfc`, `math.sqrt`, `math.exp` # and `math.log` return the correctly rounded value, so a result agrees with # CPython wherever CPython's own answer is right. import math # ── Studies as facts ───────────────────────────────────────────────────────── class Stratum: # One stratum of the study named `study`: `a_yes` successes in `a_n` # trials for arm A and `b_yes` in `b_n` for arm B. `stratum()` writes it, # named `_`. Computed: each arm's rate, and `a_wins` or # `b_wins` where that arm's rate is strictly the higher. study: string a_yes: int a_n: int b_yes: int b_n: int a_rate: float = a_yes / a_n b_rate: float = b_yes / b_n a_wins: bool = a_rate > b_rate b_wins: bool = b_rate > a_rate class Study: # A study, by a name you choose, folded from its strata: the rules sum # the successes and trials of each arm into `a_yes`, `a_n`, `b_yes` and # `b_n`. Computed from those: the pooled rates, and `pooled_a_ahead` or # `pooled_b_ahead`. The rule `paradox` concludes True where one arm is # ahead in the pool and the other arm is ahead in every stratum. The rules # also conclude, for each of the four counts, how many strata hold it # (`a_yes_strata` and so on) and how many claims of it they hold # (`a_yes_claims`), and the count is disputed where the claims are more # than the strata (`a_yes_disputed`): a stratum bound to the study twice by # that count. A disputed count is not pooled. label: string a_yes_strata: int a_n_strata: int b_yes_strata: int b_n_strata: int a_yes_claims: int a_n_claims: int b_yes_claims: int b_n_claims: int a_yes_disputed: bool a_n_disputed: bool b_yes_disputed: bool b_n_disputed: bool a_yes: int a_n: int b_yes: int b_n: int a_rate: float = a_yes / a_n b_rate: float = b_yes / b_n pooled_a_ahead: bool = a_rate > b_rate pooled_b_ahead: bool = b_rate > a_rate paradox: bool @rule def study_a_yes_strata(): for t in Study, s in Stratum: if s.study == t and s.a_yes > -1000000000000000000: t.a_yes_strata = distinct(s) @rule def study_a_n_strata(): for t in Study, s in Stratum: if s.study == t and s.a_n > -1000000000000000000: t.a_n_strata = distinct(s) @rule def study_b_yes_strata(): for t in Study, s in Stratum: if s.study == t and s.b_yes > -1000000000000000000: t.b_yes_strata = distinct(s) @rule def study_b_n_strata(): for t in Study, s in Stratum: if s.study == t and s.b_n > -1000000000000000000: t.b_n_strata = distinct(s) @rule def study_a_yes_claims(): for t in Study, s in Stratum: if s.study == t and s.a_yes > -1000000000000000000: t.a_yes_claims = count(s) @rule def study_a_n_claims(): for t in Study, s in Stratum: if s.study == t and s.a_n > -1000000000000000000: t.a_n_claims = count(s) @rule def study_b_yes_claims(): for t in Study, s in Stratum: if s.study == t and s.b_yes > -1000000000000000000: t.b_yes_claims = count(s) @rule def study_b_n_claims(): for t in Study, s in Stratum: if s.study == t and s.b_n > -1000000000000000000: t.b_n_claims = count(s) @rule def study_a_yes_disputed(): # A sum adds every value it reads, so a stratum bound to the study twice # would be counted twice, which no source gives. for t in Study: if t.a_yes_claims > t.a_yes_strata: t.a_yes_disputed = True @rule def study_a_n_disputed(): for t in Study: if t.a_n_claims > t.a_n_strata: t.a_n_disputed = True @rule def study_b_yes_disputed(): for t in Study: if t.b_yes_claims > t.b_yes_strata: t.b_yes_disputed = True @rule def study_b_n_disputed(): for t in Study: if t.b_n_claims > t.b_n_strata: t.b_n_disputed = True @rule def fold_a_yes(): for t in Study, s in Stratum: if s.study == t and absent(t.a_yes_disputed): t.a_yes = sum(s.a_yes) @rule def fold_a_n(): for t in Study, s in Stratum: if s.study == t and absent(t.a_n_disputed): t.a_n = sum(s.a_n) @rule def fold_b_yes(): for t in Study, s in Stratum: if s.study == t and absent(t.b_yes_disputed): t.b_yes = sum(s.b_yes) @rule def fold_b_n(): for t in Study, s in Stratum: if s.study == t and absent(t.b_n_disputed): t.b_n = sum(s.b_n) @rule def paradox_b_ahead_in_the_pool(): for t in Study: if t.pooled_b_ahead == True and none(s in Stratum: s.study == t and s.a_wins != True): t.paradox = True @rule def paradox_a_ahead_in_the_pool(): for t in Study: if t.pooled_a_ahead == True and none(s in Stratum: s.study == t and s.b_wins != True): t.paradox = True # ── Screening tests as facts ───────────────────────────────────────────────── class Screen: # A screening test, by a name you choose: `prevalence`, how common the # condition is, `sensitivity`, the chance a positive result is given to # someone with it, and `specificity`, the chance a negative result is # given to someone without it, each from 0 to 1. Computed by Bayes' rule: # `ppv`, the chance that a positive result is right, and `npv`, the # chance that a negative one is. `mostly_false` is True where a positive # result is more likely wrong than right. A prevalence of 0 with a # specificity of 1 has no positive results to value, and divides by zero. prevalence: float sensitivity: float specificity: float ppv: float = prevalence * sensitivity / (prevalence * sensitivity + (1.0 - prevalence) * (1.0 - specificity)) npv: float = (1.0 - prevalence) * specificity / ((1.0 - prevalence) * specificity + prevalence * (1.0 - sensitivity)) mostly_false: bool = ppv < 0.5 # ── Significance tests as facts ────────────────────────────────────────────── class Significance: # A test result, by a name you choose: its p-value and `alpha`, the level # it is held to. `significant` is True where p is below alpha. p: float alpha: float significant: bool = p < alpha # ── Writing facts ──────────────────────────────────────────────────────────── def stratum(study: string, name: string, a_yes: int, a_n: int, b_yes: int, b_n: int): # A stratum of a study, and the study itself if it is new. s = study + "_" + name update(Study(study).label, study, "strata") update(Stratum(s).study, study, "strata") update(Stratum(s).a_yes, a_yes, "strata") update(Stratum(s).a_n, a_n, "strata") update(Stratum(s).b_yes, b_yes, "strata") update(Stratum(s).b_n, b_n, "strata") def screen(name: string, prevalence: float, sensitivity: float, specificity: float): update(Screen(name).prevalence, prevalence, "screens") update(Screen(name).sensitivity, sensitivity, "screens") update(Screen(name).specificity, specificity, "screens") def significance(name: string, p: float, alpha: float): update(Significance(name).p, p, "tests") update(Significance(name).alpha, alpha, "tests") # ── Reading what the rules concluded ───────────────────────────────────────── # # Each reader reads every value in force, so a fact two sources disagree # about is reported with both values rather than refused; only a reader of a # disputed sum refuses (above). def paradoxes() -> list: # The studies whose pooled comparison reverses what every stratum shows, # sorted. out = [] vs = list_values(Study.paradox) for v in vs: parts = v.path.split(".") x = parts[len(parts) - 2] if x not in out: out = out + [x] return sorted(out) def disputed_studies() -> list: # The studies a stratum is bound to twice by one of its counts, which is # not pooled until one of its values is forgotten, sorted. out = [] for f in ["a_yes", "a_n", "b_yes", "b_n"]: if f == "a_yes": vs = list_values(Study.a_yes_disputed) elif f == "a_n": vs = list_values(Study.a_n_disputed) elif f == "b_yes": vs = list_values(Study.b_yes_disputed) else: vs = list_values(Study.b_n_disputed) for v in vs: parts = v.path.split(".") x = parts[len(parts) - 2] if x not in out: out = out + [x] return sorted(out) def pooled(study: string) -> list: # [A's pooled rate, B's pooled rate] for a study. It refuses where either # rate is absent, a study with no stratum or none with trials, and a # disputed study, which pools nothing until the dispute is settled. if study in disputed_studies(): raise ValueError("study " + study + " is disputed: a stratum of it is bound to it twice; forget the claim that is not to stand") a = current(Study.study.a_rate) b = current(Study.study.b_rate) match a.value: case real as x: match b.value: case real as y: return [x, y] case default: pass case default: pass return [] def predictive_values(name: string) -> list: # [ppv, npv] for a screening test. It refuses where either is absent. a = current(Screen.name.ppv) b = current(Screen.name.npv) match a.value: case real as x: match b.value: case real as y: return [x, y] case default: pass case default: pass return [] def significant_tests() -> list: # The significance tests whose p-value is below their level, sorted. out = [] vs = list_values(Significance.significant) for v in vs: parts = v.path.split(".") x = parts[len(parts) - 2] if v.value == True and x not in out: out = out + [x] return sorted(out) def unsettled() -> list: # Each stratum count, screen figure and test figure that holds two # different values for a fact you write, as "stratum x a_yes" and so on, # sorted. out = [] for f in ["a_yes", "a_n", "b_yes", "b_n", "prevalence", "sensitivity", "specificity", "p", "alpha"]: if f == "a_yes": vs = list_values(Stratum.a_yes) label = "stratum " elif f == "a_n": vs = list_values(Stratum.a_n) label = "stratum " elif f == "b_yes": vs = list_values(Stratum.b_yes) label = "stratum " elif f == "b_n": vs = list_values(Stratum.b_n) label = "stratum " elif f == "prevalence": vs = list_values(Screen.prevalence) label = "screen " elif f == "sensitivity": vs = list_values(Screen.sensitivity) label = "screen " elif f == "specificity": vs = list_values(Screen.specificity) label = "screen " elif f == "p": vs = list_values(Significance.p) label = "test " else: vs = list_values(Significance.alpha) label = "test " seen = {} for v in vs: parts = v.path.split(".") s = label + parts[len(parts) - 2] + " " + parts[len(parts) - 1] x = str(v.value) if s in seen and seen[s] != x and s not in out: out = out + [s] seen[s] = x return sorted(out) # ── Calculation: describing data ───────────────────────────────────────────── # # A sum or a square can pass the largest real while the statistic it serves # does not: the mean of 1e308 and 1e307, or the standard deviation of values # near 1e200. Such work is done on the values scaled by a power of two that # brings them near 1, which changes no digit of a result that does not leave # the range, and the statistic is scaled back. def scale_of(xs: list) -> int: # The power of two that brings the largest finite value of xs between 1/2 # and 1, where that value is 2^400 or more in size or under 2^-400; 0 # otherwise, where no square of the values, nor a sum of a run's worth of # them, can leave the range of a real. big = 0.0 for x in xs: if math.isfinite(x) and abs(x) > big: big = abs(x) if big == 0.0 or (big >= math.ldexp(1.0, -400) and big < math.ldexp(1.0, 400)): return 0 return math.frexp(big)[1] def mean(xs: list) -> float: # The arithmetic mean, summed left to right. xs is not empty. Where that # sum passes the largest real while every value is finite, the values are # summed again at 2^-k, k the bit length of their count, so no partial sum # can pass it. s = 0.0 for x in xs: s = s + x if math.isfinite(s): return s / len(xs) for x in xs: if not math.isfinite(x): return s / len(xs) k = math.frexp(float(len(xs)))[1] t = 0.0 for x in xs: t = t + math.ldexp(x, -k) return math.ldexp(t / len(xs), k) def variance(xs: list) -> float: # The sample variance, dividing by n - 1, by Welford's single pass over # the values scaled by scale_of(xs), and scaled back. xs holds at least # two values. e = scale_of(xs) n = 0 m = 0.0 m2 = 0.0 for v in xs: x = math.ldexp(v, -e) n = n + 1 d = x - m m = m + d / n m2 = m2 + d * (x - m) return math.ldexp(m2 / (n - 1), 2 * e) def std_dev(xs: list) -> float: # The sample standard deviation: the square root of the variance of the # values scaled by scale_of(xs), scaled back, so it is had wherever it is # a real, even where the variance is not. e = scale_of(xs) return math.ldexp(math.sqrt(variance([math.ldexp(x, -e) for x in xs])), e) def median(xs: list) -> float: # The middle value of the sorted list, or the mean of the two middle ones, # taken as the sum of their halves where their sum passes the largest real. s = sorted(xs) n = len(s) if n % 2 == 1: return s[n // 2] a = s[n // 2 - 1] b = s[n // 2] m = (a + b) / 2.0 if math.isinf(m) and math.isfinite(a) and math.isfinite(b): return a / 2.0 + b / 2.0 return m def regression(xs: list, ys: list) -> list: # [intercept, slope, r]: the least-squares line y = intercept + slope·x # and Pearson's correlation. xs has at least two distinct values, and ys # as many values as xs. The sums of products are taken of xs and ys each # scaled by its scale_of, which r does not change, and the line is scaled # back. ex = scale_of(xs) ey = scale_of(ys) sx = [math.ldexp(x, -ex) for x in xs] sy = [math.ldexp(y, -ey) for y in ys] mx = mean(sx) my = mean(sy) sxy = 0.0 sxx = 0.0 syy = 0.0 for i in range(len(sx)): sxy = sxy + (sx[i] - mx) * (sy[i] - my) sxx = sxx + (sx[i] - mx) * (sx[i] - mx) syy = syy + (sy[i] - my) * (sy[i] - my) b = sxy / sxx a = my - b * mx r = sxy / math.sqrt(sxx * syy) return [math.ldexp(a, ey), math.ldexp(b, ey - ex), r] def welch_t(xs: list, ys: list) -> float: # Welch's t statistic for two samples of unequal variance, of both samples # scaled by the one power of two scale_of finds for them together, which t # does not change. e = scale_of(xs + ys) sx = [math.ldexp(x, -e) for x in xs] sy = [math.ldexp(y, -e) for y in ys] return (mean(sx) - mean(sy)) / math.sqrt(variance(sx) / len(sx) + variance(sy) / len(sy)) def z_score(x: float, mu: float, sigma: float) -> float: # (x - mu) / sigma, taken of halves where x - mu passes the largest real. d = x - mu if math.isinf(d) and math.isfinite(x) and math.isfinite(mu): return (0.5 * x - 0.5 * mu) / sigma * 2.0 return d / sigma # ── Calculation: distributions ─────────────────────────────────────────────── def normal_pdf(x: float, mu: float, sigma: float) -> float: # The normal density. z = (x - mu) / sigma return math.exp(-0.5 * z * z) / (sigma * math.sqrt(2.0 * math.pi)) def normal_cdf(z: float) -> float: # The standard normal distribution function, 0.5·erfc(−z/√2): correctly # rounded `math.erfc`, so the tails keep their digits (z = −5 gives # 2.8665157187919455e-07). return 0.5 * math.erfc(-z / math.sqrt(2.0)) def normal_quantile(p: float) -> float: # The z whose standard normal distribution function is p, for p strictly # between 0 and 1, by 200 bisections of [−40, 40]: it is exact to the # last digit of the distribution function. lo = -40.0 hi = 40.0 for step in range(200): mid = (lo + hi) / 2.0 if normal_cdf(mid) < p: lo = mid else: hi = mid return (lo + hi) / 2.0 # The binomial and Poisson probabilities are their products wherever every # factor and step is a normal real, which is the more accurate form. Beyond, # they are Loader's saddle-point forms (C. Loader, "Fast and Accurate # Computation of Binomial Probabilities", 2000, the forms R's dbinom and dpois # take): e to a sum of terms none of which grows with n or the rate, so every # probability that is a real is had, its error growing with the size of its # logarithm. def stirling_error(n: int) -> float: # ln(n!) − ln(√(2πn)·(n/e)ⁿ), the error of Stirling's formula, for n from # 0: a table to 15 and its asymptotic series above, so no two large # logarithms are subtracted. 0 at n = 0, by Loader's convention. if n < 0: raise ValueError("n is a count, from 0") if n <= 15: t = [0.0, 0.08106146679532726, 0.0413406959554093, 0.02767792568499834, 0.020790672103765093, 0.016644691189821193, 0.013876128823070748, 0.01189670994589177, 0.010411265261972096, 0.009255462182712733, 0.00833056343336287, 0.007573675487951841, 0.00694284010720953, 0.006408994188004207, 0.0059513701127588475, 0.005554733551962801] return t[n] x = float(n) xx = x * x return (1.0 / 12.0 - (1.0 / 360.0 - (1.0 / 1260.0 - (1.0 / 1680.0 - 1.0 / 1188.0 / xx) / xx) / xx) / xx) / x def deviance(x: float, mean: float) -> float: # x·ln(x/mean) + mean − x, half the Poisson deviance of x from mean, for a # finite x from 0 and a finite mean above 0: by its series where x is # within a factor of 3 of mean, where the plain form cancels, so every # digit is kept. if not (x >= 0.0 and mean > 0.0) or math.isinf(x) or math.isinf(mean): raise ValueError("a deviance needs a finite x from 0 and a finite mean above 0") if x == 0.0: return mean d = x - mean t = x + mean if math.isinf(t): d = 0.5 * x - 0.5 * mean t = 0.5 * x + 0.5 * mean if abs(d) < 0.5 * t: v = d / t s = d * v e = 2.0 * (x * v) v = v * v for j in range(1, 200): e = e * v s1 = s + e / (2 * j + 1) if s1 == s: return s1 s = s1 return s r = x / mean if r < 2.2250738585072014e-308 or math.isinf(r): lr = math.log(x) - math.log(mean) else: lr = math.log(r) return x * lr + mean - x def exact_product(a: float, b: float) -> list: # [a·b rounded, and the part the rounding lost], which add to a·b exactly: # Dekker's product, for factors and a product under 1e300 in size. p = a * b if not (abs(a) < 1e300 and abs(b) < 1e300 and abs(p) < 1e300): raise ValueError("an exact product needs factors and a product under 1e300 in size") c = 134217729.0 * a ah = c - (c - a) al = a - ah c = 134217729.0 * b bh = c - (c - b) bl = b - bh return [p, ((ah * bh - p) + ah * bl + al * bh) + al * bl] def binom_pmf(n: int, k: int, p: float) -> float: # The probability of k successes in n trials at a chance p each, for n # from 0 and p from 0 to 1; a k outside 0 to n has probability 0. It is # C(n, k)·p^k·(1 − p)^(n − k), the coefficient a running product over the # smaller of k and n − k, and the power of 1 − p corrected by the part of # 1 − p its rounding lost, (1 − q) − p for q the rounded 1 − p, which is # exact; wherever both powers, the coefficient and the product are normal # reals. A coefficient of 515 steps or more is past the largest real, # C(1030, 515) being the first, and is not run. Otherwise it is Loader's # form, whose means n·p and n·(1 − p) are carried with the part their # rounding lost, a deviance moving by 1 − x/mean for each unit of its mean. if n < 0: raise ValueError("n is a count of trials, from 0") if not (p >= 0.0 and p <= 1.0): raise ValueError("p is a chance, from 0 to 1") if k < 0 or k > n: return 0.0 if p == 0.0: return 1.0 if k == 0 else 0.0 if p == 1.0: return 1.0 if k == n else 0.0 q = 1.0 - p fails = n - k a = p ** k b = q ** fails if a >= 2.2250738585072014e-308 and b >= 2.2250738585072014e-308 and min(k, fails) < 515: c = 1.0 for i in range(min(k, fails)): c = c * (n - i) / (i + 1) v = c * a * b * math.exp(fails * math.log1p(((1.0 - q) - p) / q)) if v >= 2.2250738585072014e-308 and not math.isinf(v): return v nf = float(n) if k == 0: return math.exp(nf * math.log1p(-p)) if k == n: return math.exp(nf * math.log(p)) kf = float(k) xf = float(n - k) m = exact_product(nf, p) r = exact_product(nf, q) r_lost = r[1] + nf * ((1.0 - q) - p) lc = (stirling_error(n) - stirling_error(k) - stirling_error(n - k) - deviance(kf, m[0]) - (1.0 - kf / m[0]) * m[1] - deviance(xf, r[0]) - (1.0 - xf / r[0]) * r_lost) return math.exp(lc - 0.5 * (1.8378770664093456 + math.log(kf) + math.log(xf / nf))) def binom_cdf(n: int, k: int, p: float) -> float: # The probability of at most k successes in n trials, summed from none. if n < 0: raise ValueError("n is a count of trials, from 0") if not (p >= 0.0 and p <= 1.0): raise ValueError("p is a chance, from 0 to 1") s = 0.0 for i in range(min(k, n) + 1): s = s + binom_pmf(n, i, p) return s def poisson_pmf(k: int, lam: float) -> float: # The probability of k events at a finite mean rate lam, from 0; a k below # 0 has probability 0. It is the running product e^−lam·lam/1·…·lam/k, # whose steps are the probabilities of 0 to k events, wherever e^−lam and # the probability are normal reals, which is for a rate to about 708; a # step falls below the smallest normal real only past the most likely # count, so the run stops there. Otherwise it is Loader's form, so a rate # past 745, where e^−lam is below the smallest real, still has the # probabilities of the counts near it. if not lam >= 0.0 or math.isinf(lam): raise ValueError("lam is a mean rate, a finite real from 0") if k < 0: return 0.0 if lam == 0.0: return 1.0 if k == 0 else 0.0 t = math.exp(-lam) if t >= 2.2250738585072014e-308: for i in range(1, k + 1): t = t * lam / i if t < 2.2250738585072014e-308: break if t >= 2.2250738585072014e-308: return t if k == 0: return math.exp(-lam) kf = float(k) return math.exp(-stirling_error(k) - deviance(kf, lam) - 0.5 * (1.8378770664093456 + math.log(kf))) def birthday(n: int) -> float: # The chance that two of n people share a birthday, 365 equally likely # days: 23 people give 0.5072972343239855. q = 1.0 for i in range(n): q = q * (365 - i) / 365.0 return 1.0 - q def bayes(prior: float, likelihood: float, false_likelihood: float) -> float: # The posterior chance of a hypothesis with the given prior, from 0 to 1, # after evidence that has `likelihood` under it and `false_likelihood` # under its negation, each a finite real from 0. Each factor is taken as # its mantissa and its power of two, so a product too small or too large # for a real still gives the posterior, which is the same to the last # digit wherever the plain products are normal reals. if not (prior >= 0.0 and prior <= 1.0): raise ValueError("the prior is a chance, from 0 to 1") if not (likelihood >= 0.0 and false_likelihood >= 0.0) or math.isinf(likelihood) or math.isinf(false_likelihood): raise ValueError("a likelihood is a finite real from 0") p = math.frexp(prior) q = math.frexp(1.0 - prior) lk = math.frexp(likelihood) fl = math.frexp(false_likelihood) a = p[0] * lk[0] b = q[0] * fl[0] return a / (a + math.ldexp(b, q[1] + fl[1] - p[1] - lk[1])) # ── Calculation: tests ─────────────────────────────────────────────────────── def chi_square(obs: list, exp: list) -> float: # Pearson's statistic, the sum of (observed − expected)² / expected, for # observed counts from 0 and expected counts above 0. A term whose square # passes the largest real, or falls below the smallest normal one, is # taken as d·(d / expected). s = 0.0 for i in range(len(obs)): if not (obs[i] >= 0.0 and exp[i] > 0.0): raise ValueError("an observed count is from 0 and an expected count above 0") d = obs[i] - exp[i] q = d * d if (math.isinf(q) or q < 2.2250738585072014e-308) and d != 0.0: s = s + d * (d / exp[i]) else: s = s + q / exp[i] return s def chi_square_p_df1(x: float) -> float: # The chance of a chi-square statistic at least x with one degree of # freedom, erfc(√(x/2)): 0.05 at 3.841458820694124. return math.erfc(math.sqrt(x / 2.0)) def chi_square_p_df2(x: float) -> float: # The same with two degrees of freedom, exp(−x/2). return math.exp(-x / 2.0) # ── Calculation: a seeded generator ────────────────────────────────────────── def lehmer_next(seed: int) -> int: # The next state of the Lehmer generator x → 48271·x mod (2³¹ − 1), for a # seed from 1 to 2³¹ − 2. The product stays under 2⁶³, so it never # overflows. A seed of 0 stays 0. return seed * 48271 % 2147483647 def lehmer_uniform(seed: int) -> float: # The state as a number between 0 and 1, the state over 2³¹ − 1. return seed / 2147483647.0 def monte_carlo_pi(n: int, seed: int) -> float: # Four times the share of n points, drawn uniformly in the unit square # from the seed, that fall inside the quarter circle. It repeats for the # same n and seed; 10,000 points from seed 1 give 3.1344. The points count # against the memory one run may hold, so a run on the service refuses at # 20,000. s = seed inside = 0 for i in range(n): s = lehmer_next(s) u = lehmer_uniform(s) s = lehmer_next(s) v = lehmer_uniform(s) if u * u + v * v < 1.0: inside = inside + 1 return 4.0 * inside / n