""" zipf-in-books.py Science Journaling Club, Volume 2 Issue 4, Summer 2026. QUESTION -------- Zipf's law says that if you rank the words of a text by how often they occur, the frequency of the word at rank r falls off as r**(-a) with a close to 1: the commonest word appears about twice as often as the second, three times as often as the third. Does that hold across languages, authors and centuries? Does the exponent a mean anything, in the sense of being stable enough that a difference between two texts is a fact about the texts rather than a fact about sampling? And where does the law fail? WHAT THIS PROGRAM IS -------------------- The club has no laboratory. The experiment here is arithmetic performed on public-domain books. Nothing was observed in a library, a classroom, or any corpus of speech. Every number below is produced by this file from fifty plain-text files downloaded from Project Gutenberg and cached under analysis/data/gutenberg. Where the text says "measured" it means "counted from those files by this code". MODEL ----- Two families are fitted, both by maximum likelihood, never by regression on a log-log plot. (1) Rank side. A text of N tokens drawn from V word types is modelled as N independent draws from a categorical distribution whose rank-r probability is p(r) = (r + q)**(-a) / Z(a, q), Z(a, q) = sum_{r=1..V} (r + q)**(-a) q = 0 is the pure Zipf law; q free is the Zipf-Mandelbrot form. Both are fitted by maximising the exact multinomial log-likelihood and compared by a likelihood-ratio test and by AIC on the same data. The ranks are estimated from the same counts that are being fitted, which biases the fit; the size and direction of that bias is measured directly in VALIDATION E rather than assumed. It turns out to push a downward, not upward as intuition about "sorting noise into an extreme order" might suggest, and the effect grows sharply as the token count shrinks. (2) Frequency side. The count x of each word type is modelled as a draw from a discrete power law p(x) = x**(-b) / zeta(b, xmin) for x >= xmin, fitted by the method of Clauset, Shalizi and Newman (2009): maximum likelihood for b at every candidate xmin, xmin chosen to minimise the Kolmogorov-Smirnov distance, and a goodness-of-fit p-value from a parametric bootstrap in which xmin is re-selected on every synthetic sample. This is the side on which an honest goodness-of-fit test is available. The two exponents are related by a = 1 / (b - 1) when both descriptions hold. ASSUMPTIONS ----------- * A token is a maximal run of Unicode letters, optionally containing internal apostrophes, lowercased. Numerals, punctuation and markup are discarded. Section 11 of the output measures what two other tokenisers do to the answer. * Word forms, not lemmas. "walk", "walks" and "walked" are three types. This matters enormously for morphologically rich languages and is the single largest modelling choice in the study. * Tokens are treated as independent draws. They are not: real text is correlated at every scale. The likelihood is therefore a working approximation and the formal standard errors it returns are too small. Section 9 measures how much too small by splitting each book in half. * Project Gutenberg boilerplate is stripped between the *** START *** and *** END *** markers. Transcriber's notes, tables of contents and translator prefaces inside those markers are left in. * Translations count as texts in the language they are written in, not in the language of the original. The Esperanto and Greek entries are translations and are labelled as such. LIMITATIONS ----------- * Fifty books is fifty books. Moreno-Sanchez et al. (2016) used 31,075. * One edition per work. Editorial choices in the Gutenberg transcription (spelling modernisation, chapter headings, end matter) are not controlled. * The goodness-of-fit machinery tests the frequency-side power law only. No comparable exact test is applied to the rank-side Zipf-Mandelbrot fit; the likelihood ratio between the two rank-side models is a relative comparison and says nothing about whether either is adequate in absolute terms. * Nothing here explains Zipf's law. The study measures it and says where it breaks. The mechanism is still open. * Bug found and fixed during validation: the first version of the discrete KS statistic (dpl_ks and the xmin search inside dpl_fit) compared the empirical CDF just below each observed value against the model CDF AT that value, instead of the model CDF just below it. At x = xmin this manufactured a fake discrepancy equal to the whole first-bin probability mass, which is exactly what VALIDATION B caught (D = 0.127 against a critical value of 0.001, on data drawn from the model being tested). The fix evaluates the model CDF at x-1 for that half of the statistic. All fifty per-book fits, xmin choices and GOF p-values below are from the corrected code. DATA ---- Source: Project Gutenberg, https://www.gutenberg.org/ Per-work URL: https://www.gutenberg.org/cache/epub//pg.txt Catalogue: https://www.gutenberg.org/cache/epub/feeds/pg_catalog.csv Cache: analysis/data/gutenberg/pg.txt, one plain UTF-8 file per work, plus analysis/data/gutenberg/MANIFEST.csv recording for every work the Gutenberg ID, title, author, language, Gutenberg release date, the URL used, the retrieval date, the byte count and the SHA-256 of the bytes. Retrieval date of the cached corpus used for the published run is printed in the output as the RETRIEVED line and stored in the MANIFEST. Project Gutenberg is a volunteer digitisation project; the works themselves are in the public domain in the United States. SEED ---- Master seed 20260614, hard-coded below. Every random stream is spawned from it through numpy SeedSequence, so the whole output is deterministic. RUN --- python zipf-in-books.py > zipf-in-books-output.txt Needs Python 3.12 and numpy. No other package is imported. The first run downloads about 50 MB from gutenberg.org; later runs read the cache and do no networking at all. """ import hashlib import math import os import re import sys import time import urllib.request from collections import Counter from datetime import date, datetime, timezone import numpy as np # Book titles and headers include non-ASCII text (French, German, Greek, # Hungarian ...). Force UTF-8 on stdout so the script runs the same way # whether the console codepage is UTF-8 or (on Windows, by default) cp1252. try: sys.stdout.reconfigure(encoding="utf-8") except AttributeError: pass MASTER_SEED = 20260614 HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "data", "gutenberg") MANIFEST = os.path.join(DATA, "MANIFEST.csv") UA = "Mozilla/5.0 (compatible; ScienceJournalingClub/1.0; educational research)" T0 = time.time() def rule(ch="=", n=78): print(ch * n) def head(title): print() rule() print(title) rule() # --------------------------------------------------------------------------- # 1. Hurwitz zeta by Euler-Maclaurin. numpy only, no scipy. # --------------------------------------------------------------------------- def hurwitz_zeta(s, q, M=16): """zeta(s, q) = sum_{k=0}^inf (q+k)**(-s) for s > 1, q > 0. Broadcasts.""" s, q = np.broadcast_arrays(np.asarray(s, dtype=float), np.asarray(q, dtype=float)) total = np.sum((q[..., None] + np.arange(M, dtype=float)) ** (-s[..., None]), axis=-1) a = q + M total = total + a ** (1.0 - s) / (s - 1.0) + 0.5 * a ** (-s) # Bernoulli corrections B2/2!, B4/4!, B6/6!, B8/8! with rising factorials total = total + s * a ** (-s - 1.0) / 12.0 total = total - s * (s + 1) * (s + 2) * a ** (-s - 3.0) / 720.0 total = total + s * (s + 1) * (s + 2) * (s + 3) * (s + 4) * a ** (-s - 5.0) / 30240.0 total = total - (s * (s + 1) * (s + 2) * (s + 3) * (s + 4) * (s + 5) * (s + 6) * a ** (-s - 7.0) / 1209600.0) return total # --------------------------------------------------------------------------- # 2. Discrete power law: p(x) = x**(-b) / zeta(b, xmin), x = xmin, xmin+1, ... # --------------------------------------------------------------------------- def dpl_mle(sum_log_x, n, xmin, lo=1.01, hi=6.0, iters=45): """Vectorised golden-section maximiser of the discrete power-law likelihood. sum_log_x, n and xmin may be arrays, one entry per independent sample.""" sum_log_x = np.asarray(sum_log_x, dtype=float) n = np.asarray(n, dtype=float) xmin = np.asarray(xmin, dtype=float) def negll(b): return n * np.log(hurwitz_zeta(b, xmin)) + b * sum_log_x gr = (math.sqrt(5.0) - 1.0) / 2.0 shape = np.broadcast(sum_log_x, n, xmin).shape a = np.full(shape, lo, dtype=float) d = np.full(shape, hi, dtype=float) c = d - gr * (d - a) e = a + gr * (d - a) fc, fe = negll(c), negll(e) for _ in range(iters): m = fc < fe d = np.where(m, e, d) a = np.where(m, a, c) c = d - gr * (d - a) e = a + gr * (d - a) fc, fe = negll(c), negll(e) return 0.5 * (a + d) def dpl_se(b, n, xmin, h=1e-4): """Asymptotic SE of b from the observed information, -d2/db2 of the loglik.""" lz = lambda x: np.log(hurwitz_zeta(x, xmin)) d2 = (lz(b + h) - 2.0 * lz(b) + lz(b - h)) / (h * h) return 1.0 / np.sqrt(np.maximum(np.asarray(n, dtype=float) * d2, 1e-300)) def dpl_cdf_at(xs, b, xmin): """Model CDF P(X <= x) for the discrete power law, evaluated at integers xs.""" return 1.0 - hurwitz_zeta(b, np.asarray(xs, dtype=float) + 1.0) / hurwitz_zeta(b, float(xmin)) def dpl_ks(tail_sorted, b, xmin): """Two-sided discrete KS distance between an ascending-sorted tail and the fit.""" n = tail_sorted.size if n < 2: return np.inf xs, idx = np.unique(tail_sorted, return_index=True) upper = np.append(idx[1:], n) / n # empirical CDF at each distinct x lower = idx / n # empirical CDF just below x cdf_at = dpl_cdf_at(xs, b, xmin) # model CDF at x, for the upper comparison cdf_below = dpl_cdf_at(xs - 1, b, xmin) # model CDF at x-1, for the lower comparison return float(max(np.max(np.abs(upper - cdf_at)), np.max(np.abs(lower - cdf_below)))) def dpl_prep(counts): """Pre-sort a sample once so the xmin sweep costs nothing extra.""" cs = np.sort(np.asarray(counts, dtype=np.int64)) xs, first = np.unique(cs, return_index=True) # first[i] = number of values < xs[i] logs = np.log(cs.astype(float)) suffix = np.concatenate([np.cumsum(logs[::-1])[::-1], [0.0]]) return cs, xs, first, suffix def dpl_fit(counts, max_cand=30, min_tail=200, prepped=None): """Clauset fit: MLE of b at every candidate xmin, xmin chosen by minimum KS. Candidates are the distinct data values, low to high, capped at max_cand and at tails of fewer than min_tail points. The floor on the tail size matters: the KS distance shrinks like 1/sqrt(ntail) all on its own, so an unbounded search walks xmin up until the tail is too small to contradict anything. The MLE runs vectorised across all candidates at once.""" cs, xs, first, suffix = dpl_prep(counts) if prepped is None else prepped n = cs.size cand = np.nonzero((n - first) >= min_tail)[0][:max_cand] if cand.size == 0: return None starts = first[cand] ntail = (n - starts).astype(float) b = dpl_mle(suffix[starts], ntail, xs[cand].astype(float)) bestD, bestk = np.inf, -1 for k, j in enumerate(cand): s = first[j] nt = float(n - s) idx = (first[j:] - s).astype(float) upper = np.append(idx[1:], nt) / nt lower = idx / nt cdf_at = dpl_cdf_at(xs[j:], b[k], float(xs[j])) cdf_below = dpl_cdf_at(xs[j:] - 1, b[k], float(xs[j])) D = max(np.max(np.abs(upper - cdf_at)), np.max(np.abs(lower - cdf_below))) if D < bestD: bestD, bestk = D, k return (int(xs[cand[bestk]]), float(b[bestk]), float(bestD), int(ntail[bestk])) class DPLSampler: """Draws from p(x) = x**(-b)/zeta(b, xmin): exact CDF inversion inside a table of K values above xmin, continuous approximation in the far tail beyond it.""" def __init__(self, b, xmin, K=6000): self.b, self.xmin, self.K = float(b), int(xmin), int(K) self.xs = np.arange(xmin, xmin + K, dtype=float) self.surv = hurwitz_zeta(self.b, self.xs) / hurwitz_zeta(self.b, float(xmin)) def draw(self, rng, n): """Inverse transform on the survival function S(x) = P(X >= x). A draw S lands on the value x with S(x) >= S > S(x+1).""" if n <= 0: return np.empty(0, dtype=np.int64) S = 1.0 - rng.random(n) # uniform on (0, 1] j = np.searchsorted(self.surv[::-1], S, side="left") far = j == 0 # beyond the end of the table out = np.empty(n, dtype=np.int64) idx = self.K - 1 - j[~far] out[~far] = self.xs[idx].astype(np.int64) if far.any(): out[far] = np.floor((self.xmin - 0.5) * S[far] ** (-1.0 / (self.b - 1.0)) + 0.5).astype(np.int64) return out def dpl_gof(rng, counts, reps=200, fit=None): """Parametric bootstrap goodness-of-fit p-value. xmin is re-selected on every synthetic sample, so the p-value pays for the xmin search. Semi-parametric: values below xmin are resampled from the empirical body of the data.""" cs = np.sort(np.asarray(counts, dtype=np.int64)) if fit is None: fit = dpl_fit(cs) if fit is None: return None xmin, b, D, ntail = fit body = cs[cs < xmin] ntot = cs.size frac_tail = ntail / ntot sampler = DPLSampler(b, xmin) nworse = 0 for _ in range(reps): k = int(rng.binomial(ntot, frac_tail)) parts = [sampler.draw(rng, k)] if ntot - k > 0 and body.size: parts.append(rng.choice(body, size=ntot - k, replace=True)) f2 = dpl_fit(np.concatenate(parts)) if f2 is not None and f2[2] >= D: nworse += 1 return dict(xmin=xmin, b=b, D=D, ntail=ntail, p=nworse / reps, reps=reps, se=float(dpl_se(b, float(ntail), float(xmin)))) # --------------------------------------------------------------------------- # 3. Rank side: Zipf and Zipf-Mandelbrot by exact multinomial likelihood. # --------------------------------------------------------------------------- def zm_loglik(counts_desc, a, q): r = np.arange(1, counts_desc.size + 1, dtype=float) Z = np.sum((r + q) ** (-a)) return float(-a * np.sum(counts_desc * np.log(r + q)) - counts_desc.sum() * math.log(Z)) def _golden_max(f, lo, hi, iters=70): gr = (math.sqrt(5.0) - 1.0) / 2.0 a, d = float(lo), float(hi) c, e = d - gr * (d - a), a + gr * (d - a) fc, fe = f(c), f(e) for _ in range(iters): if fc > fe: d, e, fe = e, c, fc c = d - gr * (d - a) fc = f(c) else: a, c, fc = c, e, fe e = a + gr * (d - a) fe = f(e) return 0.5 * (a + d) def fit_zipf(counts_desc, lo=0.3, hi=2.5): a = _golden_max(lambda x: zm_loglik(counts_desc, x, 0.0), lo, hi) return a, zm_loglik(counts_desc, a, 0.0) def fit_zm(counts_desc, qlo=-0.9, qhi=400.0): def prof(q): a = _golden_max(lambda x: zm_loglik(counts_desc, x, q), 0.3, 3.0, iters=40) return zm_loglik(counts_desc, a, q) q = _golden_max(prof, qlo, qhi, iters=50) a = _golden_max(lambda x: zm_loglik(counts_desc, x, q), 0.3, 3.0) return a, q, zm_loglik(counts_desc, a, q) def zipf_se(counts_desc, a, h=2e-4): l0 = zm_loglik(counts_desc, a, 0.0) lp = zm_loglik(counts_desc, a + h, 0.0) lm = zm_loglik(counts_desc, a - h, 0.0) d2 = -(lp - 2.0 * l0 + lm) / (h * h) return 1.0 / math.sqrt(max(d2, 1e-300)) def fit_zipf_window(counts_desc, r1, r2): """Pure Zipf exponent on ranks r1..r2 only, conditional on lying in the window.""" r2 = min(r2, counts_desc.size) if r2 - r1 + 1 < 20: return None c = counts_desc[r1 - 1:r2].astype(float) r = np.arange(r1, r2 + 1, dtype=float) tot = float(c.sum()) slog = float(np.sum(c * np.log(r))) def ll(a): return -a * slog - tot * math.log(float(np.sum(r ** (-a)))) return _golden_max(ll, 0.2, 4.0) # --------------------------------------------------------------------------- # 4. Chi-square tail probability (Wilson-Hilferty), normal tail. numpy only. # --------------------------------------------------------------------------- def norm_sf(z): return 0.5 * math.erfc(z / math.sqrt(2.0)) def chi2_sf(x, df): """Upper tail of chi-square by the Wilson-Hilferty cube-root transform.""" if df <= 0: return float("nan") z = ((x / df) ** (1.0 / 3.0) - (1.0 - 2.0 / (9.0 * df))) / math.sqrt(2.0 / (9.0 * df)) return norm_sf(z) # --------------------------------------------------------------------------- # 5. The corpus. Fifty works, fourteen written languages, sixteenth century to # twentieth, plus two classical Latin texts. The year is the conventional # first-publication date of the work as given by standard references; it is # an editorial attribution, not something this program measured. "T" marks a # translation, whose language-side text is later than the year shown. # Title, author, language and Gutenberg release date are read out of the # downloaded file itself, not typed in here. # --------------------------------------------------------------------------- CORPUS = [ # (Gutenberg ID, language tag, first published, translated?, short label) (1122, "en", 1601, False, "Hamlet"), (20, "en", 1667, False, "Paradise Lost"), (521, "en", 1719, False, "Robinson Crusoe"), (829, "en", 1726, False, "Gulliver's Travels"), (1342, "en", 1813, False, "Pride and Prejudice"), (158, "en", 1815, False, "Emma"), (84, "en", 1818, False, "Frankenstein"), (2701, "en", 1851, False, "Moby-Dick"), (98, "en", 1859, False, "A Tale of Two Cities"), (1228, "en", 1859, False, "Origin of Species"), (1400, "en", 1861, False, "Great Expectations"), (11, "en", 1865, False, "Alice in Wonderland"), (74, "en", 1876, False, "Tom Sawyer"), (76, "en", 1884, False, "Huckleberry Finn"), (1661, "en", 1892, False, "Sherlock Holmes"), (345, "en", 1897, False, "Dracula"), (219, "en", 1899, False, "Heart of Darkness"), (4300, "en", 1922, False, "Ulysses"), (64317, "en", 1925, False, "The Great Gatsby"), (6318, "fr", 1668, False, "L'Avare"), (4650, "fr", 1759, False, "Candide"), (13951, "fr", 1844, False, "Les trois mousquetaires"), (14155, "fr", 1857, False, "Madame Bovary"), (6099, "fr", 1857, False, "Les Fleurs du Mal"), (17489, "fr", 1862, False, "Les miserables I"), (5097, "fr", 1870, False, "Vingt mille lieues"), (2650, "fr", 1913, False, "Du cote de chez Swann"), (2407, "de", 1774, False, "Werther I"), (2229, "de", 1808, False, "Faust I"), (7205, "de", 1883, False, "Also sprach Zarathustra"), (5323, "de", 1895, False, "Effi Briest"), (22367, "de", 1915, False, "Die Verwandlung"), (320, "es", 1554, False, "Lazarillo de Tormes"), (2000, "es", 1605, False, "Don Quijote"), (17340, "es", 1878, False, "Marianela"), (1000, "it", 1320, False, "Divina Commedia"), (52484, "it", 1883, False, "Pinocchio"), (54829, "pt", 1881, False, "Bras Cubas"), (55752, "pt", 1899, False, "Dom Casmurro"), (15975, "nl", 1839, False, "Camera Obscura"), (11024, "nl", 1860, False, "Max Havelaar"), (218, "la", -50, False, "De Bello Gallico I-IV"), (227, "la", -19, False, "Aeneis"), (7000, "fi", 1849, False, "Kalevala"), (11940, "fi", 1870, False, "Seitseman veljesta"), (30078, "sv", 1887, False, "Hemsoborna"), (10686, "da", 1889, False, "Tine"), (43777, "hu", 1863, False, "Az uj foldesur I"), (17839, "el", -429, True, "Oidipous Tyrannos (mod. Greek trans.)"), (11511, "eo", 1719, True, "Robinsono Kruso (Esperanto trans.)"), ] LANGNAME = { "en": "English", "fr": "French", "de": "German", "es": "Spanish", "it": "Italian", "pt": "Portuguese", "nl": "Dutch", "la": "Latin", "fi": "Finnish", "sv": "Swedish", "da": "Danish", "hu": "Hungarian", "el": "Greek", "eo": "Esperanto", } def fetch(gid): """Return (raw_bytes, url, retrieved_iso, from_cache).""" os.makedirs(DATA, exist_ok=True) path = os.path.join(DATA, "pg%d.txt" % gid) stamp = os.path.join(DATA, "pg%d.when" % gid) url = "https://www.gutenberg.org/cache/epub/%d/pg%d.txt" % (gid, gid) if os.path.exists(path) and os.path.getsize(path) > 1000: with open(path, "rb") as fh: raw = fh.read() when = open(stamp).read().strip() if os.path.exists(stamp) else "unknown" return raw, url, when, True req = urllib.request.Request(url, headers={"User-Agent": UA}) raw = urllib.request.urlopen(req, timeout=120).read() with open(path, "wb") as fh: fh.write(raw) when = datetime.now(timezone.utc).date().isoformat() with open(stamp, "w") as fh: fh.write(when) time.sleep(1.0) # be polite to gutenberg.org return raw, url, when, False START_RE = re.compile(r"\*\*\*\s*START OF (?:THE|THIS) PROJECT GUTENBERG EBOOK.*?\*\*\*", re.I | re.S) END_RE = re.compile(r"\*\*\*\s*END OF (?:THE|THIS) PROJECT GUTENBERG EBOOK.*?\*\*\*", re.I | re.S) def strip_boilerplate(text): m = START_RE.search(text) if m: text = text[m.end():] m = END_RE.search(text) if m: text = text[:m.start()] return text def header_field(text, name): m = re.search(r"^%s:\s*(.+?)\s*$" % name, text[:6000], re.I | re.M) return m.group(1) if m else "" TOKEN_A = re.compile(r"[^\W\d_]+(?:['\u2019][^\W\d_]+)*", re.UNICODE) TOKEN_B = re.compile(r"[^\W\d_]+", re.UNICODE) TOKEN_C = re.compile(r"[^\W\d_]+(?:['\u2019\-][^\W\d_]+)*", re.UNICODE) def tokenise(text, pattern=TOKEN_A, fold=True): s = text.lower() if fold else text return pattern.findall(s) def counts_from_tokens(tokens): """Return counts sorted descending, plus the vocabulary size.""" tally = Counter(tokens) cnt = np.array(sorted(tally.values(), reverse=True), dtype=np.int64) return cnt, len(tally) ENGLISH_MARKERS = ("the", "and", "of", "to", "that", "with", "which", "this", "from") def english_share(tokens): n = len(tokens) if n == 0: return 1.0 hit = sum(1 for t in tokens if t in ENGLISH_MARKERS) return hit / n # --------------------------------------------------------------------------- # 6. VALIDATION. Nothing touches a real book until the machinery has been # shown to recover a known answer from data drawn from a known law. # --------------------------------------------------------------------------- def validation_zeta(): head("VALIDATION A. The Hurwitz zeta function against closed forms") print("The whole likelihood rests on zeta(s, q). It is computed here by") print("Euler-Maclaurin with 16 direct terms and four Bernoulli corrections.") print("Club value beside the analytic value, for cases where one exists.") print() apery = 1.2020569031595942854 zeta32 = 2.6123753486854883 cases = [ ("zeta(2,1) = pi^2/6", float(hurwitz_zeta(2.0, 1.0)), math.pi ** 2 / 6.0), ("zeta(4,1) = pi^4/90", float(hurwitz_zeta(4.0, 1.0)), math.pi ** 4 / 90.0), ("zeta(6,1) = pi^6/945", float(hurwitz_zeta(6.0, 1.0)), math.pi ** 6 / 945.0), ("zeta(3,1) = Apery", float(hurwitz_zeta(3.0, 1.0)), apery), ("zeta(1.5,1)", float(hurwitz_zeta(1.5, 1.0)), zeta32), ("zeta(2,5) = pi^2/6 - sum_{k=1..4} k^-2", float(hurwitz_zeta(2.0, 5.0)), math.pi ** 2 / 6.0 - sum(1.0 / k ** 2 for k in range(1, 5))), ] print("%-40s %-20s %-20s %s" % ("case", "club", "accepted", "difference")) worst = 0.0 for name, got, want in cases: d = got - want worst = max(worst, abs(d / want)) print("%-40s %-20.15f %-20.15f %+.3e" % (name, got, want, d)) # brute force check at a non-integer argument where no closed form exists s, q = 2.5, 7.5 k = np.arange(0, 4_000_000, dtype=float) brute = float(np.sum((q + k) ** (-s)) + (q + 4_000_000) ** (1 - s) / (s - 1)) got = float(hurwitz_zeta(s, q)) print("%-40s %-20.15f %-20.15f %+.3e" % ("zeta(2.5,7.5) vs 4e6-term sum", got, brute, got - brute)) worst = max(worst, abs((got - brute) / brute)) print() print("Worst relative error over all cases: %.2e" % worst) print("VERDICT: %s" % ("PASS, agreement to better than 1e-12" if worst < 1e-12 else "FAIL")) return worst def validation_sampler(seed): head("VALIDATION B. The discrete power-law sampler against its own CDF") rng = np.random.default_rng(seed) b, xmin, n = 1.95, 7, 2_000_000 x = np.sort(DPLSampler(b, xmin).draw(rng, n)) D = dpl_ks(x, b, xmin) crit = 1.358 / math.sqrt(n) print("b = %.2f, xmin = %d, n = %s draws" % (b, xmin, format(n, ","))) print("KS distance from the exact model CDF : %.6f" % D) print("95%% KS critical value 1.358/sqrt(n) : %.6f" % crit) print("VERDICT: %s" % ("PASS, the sampler is drawing from the law it claims" if D < crit else "FAIL")) print() b2, xmin2 = 2.5, 7 x2 = DPLSampler(b2, xmin2).draw(rng, 4_000_000) got = float(x2.mean()) want = float(hurwitz_zeta(b2 - 1.0, float(xmin2)) / hurwitz_zeta(b2, float(xmin2))) se = float(x2.std(ddof=1) / math.sqrt(x2.size)) print("Analytic mean at b = 2.5, xmin = 7 is zeta(1.5,7)/zeta(2.5,7).") print("club mean %.5f +/- %.5f, analytic %.5f, difference %+.5f = %.2f SE" % (got, se, want, got - want, (got - want) / se)) print("VERDICT: %s" % ("PASS" if abs(got - want) < 4 * se else "FAIL")) return D, crit def validation_recovery(seed): head("VALIDATION C. Recovering a known exponent from synthetic data") rng = np.random.default_rng(seed) print("Draws from a discrete power law with a known b, fitted by the same MLE") print("used on the books. xmin is held at its true value here so that the") print("estimator is tested on its own, without the xmin search on top of it.") print("reps = 300 per cell. z is the bias divided by its own standard error;") print("|z| > 3 would mean a real bias rather than Monte Carlo noise.") print() print("%-7s %-8s %-10s %-11s %-10s %-10s %s" % ("b_true", "n", "mean b", "bias", "SD(b)", "mean SE", "z")) reps = 300 conv = {} worst_z = 0.0 for b_true in (1.6, 1.95, 2.3, 3.0): sampler = DPLSampler(b_true, 7) for n in (250, 1000, 4000, 16000, 64000): sl = np.empty(reps) for i in range(reps): sl[i] = np.log(sampler.draw(rng, n).astype(float)).sum() bh = dpl_mle(sl, np.full(reps, float(n)), np.full(reps, 7.0)) bias = float(bh.mean() - b_true) sd = float(bh.std(ddof=1)) mse = float(np.mean(dpl_se(bh, float(n), 7.0))) z = bias / (sd / math.sqrt(reps)) worst_z = max(worst_z, abs(z)) print("%-7.2f %-8d %-10.4f %+-11.4f %-10.4f %-10.4f %+.2f" % (b_true, n, bh.mean(), bias, sd, mse, z)) if abs(b_true - 1.95) < 1e-9: conv[n] = (float(bh.mean()), sd, mse, bias) print() print("Largest |z| over the twenty cells: %.2f" % worst_z) print("The asymptotic standard error tracks the observed scatter, which is") print("what licenses quoting SEs for the books without bootstrapping each one.") print() print("FIGDATA convergence n mean sd se") for n in sorted(conv): m, sd, mse, bias = conv[n] print("FIGDATA conv %d %.5f %.5f %.5f" % (n, m, sd, mse)) return worst_z, conv def validation_gof(seed): head("VALIDATION D. Is the goodness-of-fit test calibrated?") rng = np.random.default_rng(seed) print("A p-value is only worth printing if it is uniform when the null is true.") print("Here: 80 synthetic vocabularies of 6,000 word types whose counts really") print("are a discrete power law (b = 1.95, xmin = 1), each put through the whole") print("procedure, xmin search included, with 80 bootstrap replicates.") ntypes, b_true, reps_in, nsets = 6000, 1.95, 80, 80 sampler = DPLSampler(b_true, 1) ps = [] for _ in range(nsets): counts = sampler.draw(rng, ntypes) g = dpl_gof(rng, counts, reps=reps_in) if g: ps.append(g["p"]) ps = np.array(ps) r05 = float(np.mean(ps < 0.05)) r10 = float(np.mean(ps <= 0.10)) se05 = math.sqrt(0.05 * 0.95 / ps.size) print() print("rejection rate at p < 0.05 : %.3f (nominal 0.050, SE %.3f, z = %+.2f)" % (r05, se05, (r05 - 0.05) / se05)) print("rejection rate at p <= 0.10: %.3f (nominal 0.100)" % r10) print("mean p = %.3f (uniform would give 0.500), median p = %.3f" % (ps.mean(), np.median(ps))) print("VERDICT: %s" % ("PASS, the test does not reject data that obey the law" if abs(r05 - 0.05) < 3 * se05 else "MARGINAL, see the printed z")) print() print("And the other half of calibration: does it reject when it should?") print("40 synthetic vocabularies of the same size drawn from a log-normal,") print("which is the classic look-alike for a power law.") rej = 0 tot = 0 for _ in range(40): counts = np.maximum(1, np.round(rng.lognormal(mean=0.9, sigma=1.6, size=ntypes))).astype(np.int64) g = dpl_gof(rng, counts, reps=reps_in) if g: tot += 1 rej += (g["p"] < 0.05) print("power against log-normal counts at p < 0.05: %d of %d = %.2f" % (rej, tot, rej / max(tot, 1))) return r05, ps.size def validation_rank(seed): head("VALIDATION E. The rank-side fit, and the price of ranking your own data") rng = np.random.default_rng(seed) print("A synthetic book: N tokens drawn from a Zipf-Mandelbrot distribution over") print("V types with known a and q. The fit then does what it does to a real") print("book, which is to rank the types by their observed counts. Those ranks") print("are estimates, and the noise in them pushes the exponent up. This is the") print("measurement of that bias, not an argument that it is small.") print() out = {} for (a_true, q_true, V, N) in ((1.00, 0.0, 20000, 200000), (1.05, 2.0, 20000, 200000), (1.00, 0.0, 20000, 20000)): r = np.arange(1, V + 1, dtype=float) p = (r + q_true) ** (-a_true) p = p / p.sum() ah, qh, azh = [], [], [] reps = 20 for _ in range(reps): c = rng.multinomial(N, p) c = np.sort(c[c > 0])[::-1].astype(np.int64) a0, _ = fit_zipf(c) a1, q1, _ = fit_zm(c) ah.append(a1) qh.append(q1) azh.append(a0) ah, qh, azh = np.array(ah), np.array(qh), np.array(azh) print("truth a = %.2f, q = %.1f, V = %s, N = %s tokens, %d replicates" % (a_true, q_true, format(V, ","), format(N, ","), reps)) print(" Zipf-Mandelbrot a : %.4f +/- %.4f bias %+.4f (%+.1f%%)" % (ah.mean(), ah.std(ddof=1), ah.mean() - a_true, 100 * (ah.mean() - a_true) / a_true)) print(" Zipf-Mandelbrot q : %.3f +/- %.3f truth %.3f" % (qh.mean(), qh.std(ddof=1), q_true)) print(" pure Zipf a : %.4f +/- %.4f bias %+.4f (%+.1f%%)" % (azh.mean(), azh.std(ddof=1), azh.mean() - a_true, 100 * (azh.mean() - a_true) / a_true)) print() out[(a_true, q_true, N)] = (float(ah.mean()), float(ah.std(ddof=1)), float(azh.mean()), float(azh.std(ddof=1))) print("Read the third block against the first: with ten times fewer tokens the") print("bias grows, because rank noise grows. Every exponent quoted for a real") print("book below carries a bias of this kind and this size. It is the same") print("sign and roughly the same size for every book, so comparisons between") print("books survive it far better than any single value does.") return out # --------------------------------------------------------------------------- # 7. The corpus: download or read the cache, record provenance, tokenise. # --------------------------------------------------------------------------- def load_corpus(): head("THE CORPUS. Fifty works from Project Gutenberg") print("Source https://www.gutenberg.org/ , one plain-text file per work at") print("https://www.gutenberg.org/cache/epub//pg.txt") print("Title, author, language and release date are read from the file header.") print("Bytes are counted before Gutenberg boilerplate is stripped.") print() rows = [] manifest = ["gutenberg_id,language_tag,first_published,translation,title,author," "gutenberg_language,gutenberg_release_date,bytes,sha256,url,retrieved"] for gid, lang, year, trans, label in CORPUS: raw, url, when, cached = fetch(gid) text = raw.decode("utf-8", errors="replace") title = header_field(text, "Title") or label author = (header_field(text, "Author") or header_field(text, "Compiler") or header_field(text, "Editor") or "unattributed") gl = header_field(text, "Language") rel = header_field(text, "Release date") rel = re.sub(r"\s*\[.*$", "", rel) sha = hashlib.sha256(raw).hexdigest() body = strip_boilerplate(text) rows.append(dict(gid=gid, lang=lang, year=year, trans=trans, label=label, title=title, author=author, glang=gl, release=rel, nbytes=len(raw), sha=sha, url=url, when=when, body=body)) manifest.append('%d,%s,%d,%s,"%s","%s","%s","%s",%d,%s,%s,%s' % (gid, lang, year, "yes" if trans else "no", title.replace('"', "'"), author.replace('"', "'"), gl, rel, len(raw), sha, url, when)) os.makedirs(DATA, exist_ok=True) with open(MANIFEST, "w", encoding="utf-8") as fh: fh.write("\n".join(manifest) + "\n") whens = sorted({r["when"] for r in rows}) print("RETRIEVED %s" % (", ".join(whens))) print("MANIFEST written to analysis/data/gutenberg/MANIFEST.csv") print() print("%-4s %-7s %-6s %-34s %-26s %-13s %-10s %s" % ("#", "PG ID", "lang", "title (from the file header)", "author", "PG release", "kB", "sha256[:8]")) for i, r in enumerate(rows, 1): print("%-4d %-7d %-6s %-34s %-26s %-13s %-10.1f %s" % (i, r["gid"], r["lang"], r["title"][:34], r["author"][:26], r["release"][:13], r["nbytes"] / 1024.0, r["sha"][:8])) return rows def quality_check(rows): head("DATA QUALITY. Is each file actually in the language it claims?") print("A crude but effective gate: the share of tokens that are common English") print("function words. English books should score high, everything else low.") print("A non-English file scoring high would mean English editorial matter has") print("leaked into the text and is being counted as if it were the work.") print() print("%-34s %-6s %-10s %s" % ("title", "lang", "EN share", "flag")) flagged = [] for r in rows: share = english_share(r["tokens"]) r["en_share"] = share flag = "" if r["lang"] != "en" and share > 0.02: flag = "FLAG" flagged.append(r["label"]) print("%-34s %-6s %-10.4f %s" % (r["label"][:34], r["lang"], share, flag)) print() if flagged: print("Flagged: %s" % ", ".join(flagged)) print("Those files carry English apparatus (notes, vocabulary, a translator's") print("preface) inside the Gutenberg markers. They are kept, and every") print("language-level statement below reports whether dropping them changes it.") else: print("Nothing flagged.") return flagged # --------------------------------------------------------------------------- # 8. Per-book fits. # --------------------------------------------------------------------------- WINDOWS = [(1, 100), (101, 1000), (1001, 10000), (10001, 10 ** 9)] def analyse_books(rows, seed, gof_reps=200): head("PER-BOOK FITS") print("Rank side: pure Zipf p(r) ~ r^-a and Zipf-Mandelbrot p(r) ~ (r+q)^-a, both") print("by exact multinomial maximum likelihood over all V ranks. dAIC is") print("AIC(Zipf) - AIC(Zipf-Mandelbrot); positive means the shifted form wins.") print("Frequency side: Clauset MLE for b with xmin by minimum KS, and a") print("goodness-of-fit p from %d parametric bootstrap replicates each." % gof_reps) print("p <= 0.10 is the Clauset rule for ruling the pure power law out.") print() rng = np.random.default_rng(seed) print("%-32s %-9s %-8s %-7s %-8s %-8s %-8s %-7s %-6s %-6s %-5s %s" % ("title", "tokens", "types", "hapax", "a_zipf", "SE", "a_ZM", "q", "dAIC", "b", "xmin", "p_GOF")) for r in rows: c = r["counts"] N = int(c.sum()) V = int(c.size) hap = float(np.mean(c == 1)) a_z, ll_z = fit_zipf(c) se_z = zipf_se(c, a_z) a_m, q_m, ll_m = fit_zm(c) lr = 2.0 * (ll_m - ll_z) daic = lr - 2.0 fit = dpl_fit(c) g = dpl_gof(rng, c, reps=gof_reps, fit=fit) if fit else None r.update(N=N, V=V, hapax=hap, a_z=a_z, se_z=se_z, a_m=a_m, q=q_m, ll_z=ll_z, ll_m=ll_m, lr=lr, daic=daic, b=g["b"] if g else float("nan"), b_se=g["se"] if g else float("nan"), xmin=g["xmin"] if g else -1, ntail=g["ntail"] if g else -1, ks=g["D"] if g else float("nan"), p_gof=g["p"] if g else float("nan")) r["a_from_b"] = 1.0 / (r["b"] - 1.0) r["win"] = [fit_zipf_window(c, w0, w1) for (w0, w1) in WINDOWS] print("%-32s %-9s %-8s %-7.3f %-8.4f %-8.4f %-8.4f %-7.2f %-6.0f %-6.3f %-5d %.3f" % (r["label"][:32], format(N, ","), format(V, ","), hap, a_z, se_z, a_m, q_m, daic, r["b"], r["xmin"], r["p_gof"])) return rows def weighted_stats(vals, ses): w = 1.0 / np.asarray(ses) ** 2 mean = float(np.sum(w * vals) / np.sum(w)) Q = float(np.sum(w * (np.asarray(vals) - mean) ** 2)) df = len(vals) - 1 return mean, Q, df def group_tests(rows): head("DOES THE EXPONENT MEAN ANYTHING? Heterogeneity across the fifty books") a = np.array([r["a_z"] for r in rows]) se = np.array([r["se_z"] for r in rows]) mean, Q, df = weighted_stats(a, se) print("Unweighted mean a over 50 books : %.4f, SD %.4f, range %.4f to %.4f" % (a.mean(), a.std(ddof=1), a.min(), a.max())) print("Inverse-variance weighted mean : %.4f" % mean) print("Cochran Q : %.1f on %d degrees of freedom" % (Q, df)) print("Q/df : %.1f (1.0 would mean the books agree)" % (Q / df)) print("p : %.3g" % chi2_sf(Q, df)) print() print("The formal standard errors are tiny because the likelihood pretends every") print("token is an independent draw. They are not a believable yardstick for") print("comparing books, and section 9 replaces them with a measured one.") print() print("By language (original-language texts only where marked):") print("%-13s %-4s %-9s %-9s %-9s %s" % ("language", "n", "mean a", "SD", "min", "max")) langstat = {} for lg in sorted({r["lang"] for r in rows}): vals = np.array([r["a_z"] for r in rows if r["lang"] == lg]) langstat[lg] = vals print("%-13s %-4d %-9.4f %-9s %-9.4f %.4f" % (LANGNAME[lg], vals.size, vals.mean(), ("%.4f" % vals.std(ddof=1)) if vals.size > 1 else "n/a", vals.min(), vals.max())) print() print("Same author, two books, same language: the tightest comparison available.") pairs = [("Pride and Prejudice", "Emma"), ("A Tale of Two Cities", "Great Expectations"), ("Tom Sawyer", "Huckleberry Finn"), ("Bras Cubas", "Dom Casmurro")] by = {r["label"]: r for r in rows} out_pairs = [] for x, y in pairs: rx, ry = by[x], by[y] d = rx["a_z"] - ry["a_z"] z = d / math.sqrt(rx["se_z"] ** 2 + ry["se_z"] ** 2) out_pairs.append((x, y, d, z)) print(" %-24s vs %-22s da = %+.4f z(formal) = %+8.1f" % (x, y, d, z)) print() print("Morphology. Finnish and Hungarian inflect heavily, so a given lemma") print("appears as many distinct word forms and the head of the distribution is") print("flattened. English and Chinese-style isolating languages do the opposite.") en = langstat["en"] rich = np.concatenate([langstat["fi"], langstat["hu"]]) print(" English n=%d mean a = %.4f" % (en.size, en.mean())) print(" Finnish + Hungarian n=%d mean a = %.4f" % (rich.size, rich.mean())) print(" difference %+.4f" % (en.mean() - rich.mean())) ttl = math.sqrt(en.var(ddof=1) / en.size + rich.var(ddof=1) / rich.size) print(" Welch t = %+.2f using between-book scatter as the error" % ((en.mean() - rich.mean()) / ttl)) print() print("Type-token ratio, the same story without any fitting at all:") for lg in ("en", "fi", "hu", "la", "eo"): vv = [r["V"] / r["N"] for r in rows if r["lang"] == lg] print(" %-11s mean types per token %.4f" % (LANGNAME[lg], float(np.mean(vv)))) print() orig = [r for r in rows if not r["trans"] and r["year"] > 1400] yr = np.array([r["year"] for r in orig], dtype=float) av = np.array([r["a_z"] for r in orig]) A = np.vstack([np.ones_like(yr), yr - yr.mean()]).T coef, res, _, _ = np.linalg.lstsq(A, av, rcond=None) resid = av - A @ coef s2 = float(resid @ resid) / (yr.size - 2) cov = s2 * np.linalg.inv(A.T @ A) slope, sse = coef[1], math.sqrt(cov[1, 1]) print("Century trend, %d original-language works from %d to %d:" % (yr.size, int(yr.min()), int(yr.max()))) print(" d a / d century = %+.4f +/- %.4f (t = %+.2f, p = %.3f)" % (100 * slope, 100 * sse, slope / sse, 2 * norm_sf(abs(slope / sse)))) return dict(mean=float(a.mean()), sd=float(a.std(ddof=1)), Q=Q, df=df, wmean=mean, langstat=langstat, pairs=out_pairs, slope=100 * slope, slope_se=100 * sse) # --------------------------------------------------------------------------- # 9. Where the law fails: the tail, the shifted form, and the hapax count. # --------------------------------------------------------------------------- def tail_analysis(rows): head("FAILURE MODE 1. The exponent is not one number, it drifts with rank") print("The same pure Zipf law fitted inside four rank windows, conditional on") print("the rank falling in the window. If a single power law described the") print("whole distribution these four numbers would agree.") print() names = ["r 1-100", "r 101-1000", "r 1001-10000", "r >10000"] print("%-32s %-11s %-11s %-13s %s" % ("title", names[0], names[1], names[2], names[3])) for r in rows: w = r["win"] print("%-32s %-11s %-11s %-13s %s" % (r["label"][:32], "%.3f" % w[0] if w[0] else "n/a", "%.3f" % w[1] if w[1] else "n/a", "%.3f" % w[2] if w[2] else "n/a", "%.3f" % w[3] if w[3] else "n/a")) print() cols = [] for j in range(4): v = np.array([r["win"][j] for r in rows if r["win"][j] is not None]) cols.append(v) print("%-14s n = %-4d mean %.4f SD %.4f" % (names[j], v.size, v.mean(), v.std(ddof=1))) d = np.array([r["win"][2] - r["win"][0] for r in rows if r["win"][2] is not None and r["win"][0] is not None]) print() print("Within-book change from the head window to the r 1001-10000 window:") print(" mean %+.4f, SD %.4f, positive in %d of %d books" % (d.mean(), d.std(ddof=1), int((d > 0).sum()), d.size)) print(" paired t = %+.2f" % (d.mean() / (d.std(ddof=1) / math.sqrt(d.size)))) return cols, d def mandelbrot_analysis(rows): head("FAILURE MODE 2. A pure power law against Mandelbrot's shifted form") print("Zipf-Mandelbrot adds one parameter, the shift q, and nests the pure law") print("at q = 0. Twice the log-likelihood difference is a chi-square on 1 df.") print() lr = np.array([r["lr"] for r in rows]) wins = int((np.array([r["daic"] for r in rows]) > 0).sum()) print("Zipf-Mandelbrot beats pure Zipf on AIC in %d of %d books." % (wins, len(rows))) print("Likelihood-ratio statistic: median %.0f, min %.1f, max %.0f" % (np.median(lr), lr.min(), lr.max())) print("A chi-square on 1 df needs 3.84 for p = 0.05. %d of %d books exceed it." % (int((lr > 3.84).sum()), len(rows))) q = np.array([r["q"] for r in rows]) print("Fitted shift q: median %.2f, range %.2f to %.2f" % (np.median(q), q.min(), q.max())) print() print("What the shift does to the exponent:") az = np.array([r["a_z"] for r in rows]) am = np.array([r["a_m"] for r in rows]) print(" mean a under pure Zipf %.4f" % az.mean()) print(" mean a under Zipf-Mandelbrot %.4f" % am.mean()) print(" mean change %+.4f, SD %.4f" % ((am - az).mean(), (am - az).std(ddof=1))) print() print("So the headline exponent moves by about %.0f%% depending on which of two" % (100 * abs((am - az).mean()) / az.mean())) print("defensible models you fit. Anyone quoting a Zipf exponent to three") print("decimals without saying which form they fitted is quoting noise.") return wins, float(np.median(q)), float(am.mean() - az.mean()) def gof_summary(rows): head("FAILURE MODE 3. Does the pure power law actually fit? The GOF test") p = np.array([r["p_gof"] for r in rows]) b = np.array([r["b"] for r in rows]) print("Clauset, Shalizi and Newman rule a power law out when p <= 0.10.") print(" books with p > 0.10 (power law not ruled out): %d of %d" % (int((p > 0.10).sum()), p.size)) print(" books with p >= 0.05 : %d of %d" % (int((p >= 0.05).sum()), p.size)) print(" median p = %.3f" % np.median(p)) print() print("Frequency-side exponent b over the fifty books: mean %.3f, SD %.3f" % (b.mean(), b.std(ddof=1))) ok = p > 0.10 if ok.any(): print(" restricted to books the test does not reject: mean %.3f, SD %.3f, n = %d" % (b[ok].mean(), b[ok].std(ddof=1) if ok.sum() > 1 else float("nan"), int(ok.sum()))) ab = np.array([r["a_from_b"] for r in rows]) az = np.array([r["a_z"] for r in rows]) print() print("Two roads to the same exponent. a = 1/(b-1) from the frequency side") print("against the rank-side a, per book:") print(" mean 1/(b-1) %.4f" % ab.mean()) print(" mean rank a %.4f" % az.mean()) print(" mean difference %+.4f, SD %.4f, correlation %.3f" % ((ab - az).mean(), (ab - az).std(ddof=1), float(np.corrcoef(ab, az)[0, 1]))) print("They are measuring the same thing through different windows: the rank fit") print("uses every type, the frequency fit uses only the %d%% of types above xmin." % round(100 * np.mean([r["ntail"] / r["V"] for r in rows]))) return float(np.median(p)), int((p > 0.10).sum()) def hapax_check(rows): head("FAILURE MODE 4. Words that appear exactly once") print("Under a fitted model the expected number of types seen exactly once is") print("sum_r N p_r (1-p_r)^(N-1). Books really are full of one-off words, and") print("this is where a rank-side model is asked to predict something it was not") print("fitted on.") print() print("%-32s %-10s %-12s %-12s %s" % ("title", "observed", "Zipf pred.", "ZM pred.", "ZM/obs")) rat = [] for r in rows: V, N = r["V"], r["N"] rr = np.arange(1, V + 1, dtype=float) for key, a, q in (("z", r["a_z"], 0.0), ("m", r["a_m"], r["q"])): p = (rr + q) ** (-a) p = p / p.sum() r["hap_" + key] = float(np.sum(N * p * np.exp((N - 1) * np.log1p(-p)))) obs = float((r["counts"] == 1).sum()) r["hap_obs"] = obs rat.append(r["hap_m"] / obs) print("%-32s %-10d %-12.0f %-12.0f %.3f" % (r["label"][:32], int(obs), r["hap_z"], r["hap_m"], r["hap_m"] / obs)) rat = np.array(rat) print() print("Zipf-Mandelbrot predicted / observed hapax count: mean %.3f, SD %.3f" % (rat.mean(), rat.std(ddof=1))) print("under-predicts in %d of %d books" % (int((rat < 1).sum()), rat.size)) return float(rat.mean()) # --------------------------------------------------------------------------- # 10. How big is a real difference? Split-half, and tokeniser sensitivity. # --------------------------------------------------------------------------- def split_half(rows): head("HOW BIG IS A REAL DIFFERENCE? Splitting every book in half") print("A book is not a bag of independent draws, so the likelihood's standard") print("error understates the true variability. Measure it instead: fit the") print("first half and the second half of each book separately. The scatter") print("between halves of the same book is the instrument's real repeatability.") print() print("%-32s %-9s %-9s %-9s %s" % ("title", "a first", "a second", "diff", "formal SE(half)")) diffs, formals = [], [] for r in rows: toks = r["tokens"] h = len(toks) // 2 out = [] for part in (toks[:h], toks[h:]): c, _ = counts_from_tokens(part) a, _ = fit_zipf(c) out.append((a, zipf_se(c, a))) d = out[0][0] - out[1][0] f = math.sqrt(out[0][1] ** 2 + out[1][1] ** 2) diffs.append(d) formals.append(f) r["half_diff"] = d print("%-32s %-9.4f %-9.4f %+-9.4f %.4f" % (r["label"][:32], out[0][0], out[1][0], d, f)) diffs = np.array(diffs) formals = np.array(formals) sd_obs = float(diffs.std(ddof=1)) sd_formal = float(formals.mean()) print() print("observed SD of the half-to-half difference : %.4f" % sd_obs) print("formal SE of that same difference : %.4f" % sd_formal) print("ratio : %.1f" % (sd_obs / sd_formal)) print() sigma_emp = sd_obs / math.sqrt(2.0) print("Take %.4f as the honest one-book standard error, which is the half-to-half" % sigma_emp) print("scatter divided by sqrt(2). It is %.0f times the likelihood's own figure." % (sigma_emp / np.mean([r["se_z"] for r in rows]))) a = np.array([r["a_z"] for r in rows]) Q = float(np.sum((a - a.mean()) ** 2) / sigma_emp ** 2) df = a.size - 1 print("Redo the heterogeneity test with it: Q = %.0f on %d df, Q/df = %.0f, p = %.3g" % (Q, df, Q / df, chi2_sf(Q, df))) print("The differences between books are still far larger than the noise.") print() print("Same-author pairs against the measured scatter:") by = {r["label"]: r for r in rows} for x, y in (("Pride and Prejudice", "Emma"), ("A Tale of Two Cities", "Great Expectations"), ("Tom Sawyer", "Huckleberry Finn"), ("Bras Cubas", "Dom Casmurro")): d = by[x]["a_z"] - by[y]["a_z"] print(" %-24s vs %-22s da = %+.4f z(measured) = %+.2f" % (x, y, d, d / (sigma_emp * math.sqrt(2)))) return sigma_emp, sd_obs, Q, df def tokeniser_sensitivity(rows): head("SENSITIVITY. What a different tokeniser does to the answer") print("Tokeniser A, used everywhere above: letters with internal apostrophes,") print("lowercased. B: letters only, so don't becomes don and t. C: letters with") print("internal apostrophes and hyphens, so to-morrow stays one word. D: A") print("without lowercasing, so The and the are different types.") print() print("%-32s %-9s %-9s %-9s %s" % ("title", "A", "B - A", "C - A", "D - A")) dB, dC, dD = [], [], [] for r in rows: base = r["a_z"] vals = [] for pat, fold in ((TOKEN_B, True), (TOKEN_C, True), (TOKEN_A, False)): c, _ = counts_from_tokens(tokenise(r["body"], pat, fold)) vals.append(fit_zipf(c)[0] - base) dB.append(vals[0]) dC.append(vals[1]) dD.append(vals[2]) print("%-32s %-9.4f %+-9.4f %+-9.4f %+.4f" % (r["label"][:32], base, vals[0], vals[1], vals[2])) for name, v in (("B letters only", dB), ("C hyphens joined", dC), ("D case kept", dD)): v = np.array(v) print("%-20s mean %+.4f SD %.4f largest |change| %.4f" % (name, v.mean(), v.std(ddof=1), np.abs(v).max())) print() print("Set those against the %.4f half-to-half scatter: the choice of tokeniser" % SIGMA_EMP[0]) print("moves the exponent by more than the measurement noise does, which is why") print("comparing an exponent from one study with an exponent from another is") print("only meaningful when both tokenised the same way.") return np.array(dB), np.array(dC), np.array(dD) SIGMA_EMP = [float("nan")] # --------------------------------------------------------------------------- # 11. Against the published literature. # --------------------------------------------------------------------------- def published_comparison(rows): head("CLUB VALUES BESIDE PUBLISHED VALUES") by = {r["label"]: r for r in rows} md = by["Moby-Dick"] print("1. Moby-Dick, the single most-fitted word-frequency data set in the") print(" power-law literature. Clauset, Shalizi and Newman (2009), Table 6.1,") print(" row 'count of word use', report for exactly this book:") print(" n(types) = 18,855 xmin = 7 +/- 2 b = 1.95 +/- 0.02") print(" ntail = 2958 +/- 987 p = 0.49") print(" The club, from its own download and its own tokeniser:") print(" n(types) = %s xmin = %d b = %.3f +/- %.3f" % (format(md["V"], ","), md["xmin"], md["b"], md["b_se"])) print(" ntail = %d p = %.3f KS = %.4f" % (md["ntail"], md["p_gof"], md["ks"])) d = md["b"] - 1.95 comb = math.sqrt(md["b_se"] ** 2 + 0.02 ** 2) print(" difference in b: %+.3f, which is %.2f combined standard errors." % (d, d / comb)) print(" %s" % ("AGREEMENT." if abs(d / comb) < 3 else "DISAGREEMENT at %.1f sigma. See the note below." % abs(d / comb))) print() print(" The two studies disagree about xmin, %d against 7, and that is the more" % md["xmin"]) print(" interesting disagreement, so here is the same fit held at their xmin on") print(" our own token counts, which is the closest thing to a like-for-like") print(" comparison the two studies admit:") cs = np.sort(md["counts"]) tail = cs[np.searchsorted(cs, 7):] b7 = float(dpl_mle(float(np.log(tail).sum()), float(tail.size), 7.0)) se7 = float(dpl_se(b7, float(tail.size), 7.0)) d7 = b7 - 1.95 c7 = math.sqrt(se7 ** 2 + 0.02 ** 2) print(" club at xmin = 7: ntail = %d, b = %.4f +/- %.4f" % (tail.size, b7, se7)) print(" Clauset et al. : ntail = 2958 +/- 987, b = 1.95 +/- 0.02") print(" difference %+.4f, which is %.2f combined standard errors. %s" % (d7, d7 / c7, "AGREEMENT." if abs(d7 / c7) < 3 else "DISAGREEMENT.")) print(" Our KS search lands on xmin = %d, one word past their %d; both are inside" % (md["xmin"], 7)) print(" the +/- 2 they quote for their own bootstrap, so the two studies agree") print(" that this parameter is the soft one. The exponent is not: at a common") print(" xmin the two agree to %.1f standard errors." % abs(d7 / c7)) print() print(" Note on what can and cannot be compared here. Clauset et al. do not") print(" publish their tokeniser. Our vocabulary is %s types against their" % format(md["V"], ",")) print(" 18,855, a difference of %+.1f%%, so the two studies are not counting" % (100 * (md["V"] - 18855) / 18855)) print(" quite the same objects and a small shift in b is expected on those") print(" grounds alone. The comparison is worth making anyway because it is") print(" the one number in this study that somebody else has published for the") print(" same book.") print() b = np.array([r["b"] for r in rows]) ben = np.array([r["b"] for r in rows if r["lang"] == "en"]) print("2. Moreno-Sanchez, Font-Clos and Corral (2016) fitted 31,075 English") print(" Gutenberg texts with the same MLE-plus-KS machinery and report a mean") print(" frequency-side exponent of about 2.03 with SD 0.15 among the texts") print(" their test did not reject, and roughly 40% of texts surviving the test.") print(" Club, English subset, n = %d: mean b = %.3f, SD %.3f" % (ben.size, ben.mean(), ben.std(ddof=1))) print(" Club, all %d books: mean b = %.3f, SD %.3f" % (b.size, b.mean(), b.std(ddof=1))) se_mean = ben.std(ddof=1) / math.sqrt(ben.size) print(" difference from 2.03 on the English subset: %+.3f = %.2f SE of our mean" % (ben.mean() - 2.03, (ben.mean() - 2.03) / se_mean)) p = np.array([r["p_gof"] for r in rows]) pen = np.array([r["p_gof"] for r in rows if r["lang"] == "en"]) print(" Club survival rate at p > 0.10: %d of %d English (%.0f%%), %d of %d overall (%.0f%%)" % (int((pen > 0.1).sum()), pen.size, 100 * np.mean(pen > 0.1), int((p > 0.1).sum()), p.size, 100 * np.mean(p > 0.1))) print() az = np.array([r["a_z"] for r in rows]) aen = np.array([r["a_z"] for r in rows if r["lang"] == "en"]) print("3. Zipf's own claim, and the textbook figure, is a rank exponent of about") print(" 1. Equivalently b = 1 + 1/a = 2.") print(" Club rank exponent, English books : mean %.4f, SD %.4f" % (aen.mean(), aen.std(ddof=1))) print(" Club rank exponent, all books : mean %.4f, SD %.4f" % (az.mean(), az.std(ddof=1))) print(" distance of the English mean from 1: %+.4f" % (aen.mean() - 1.0)) print(" in units of the measured one-book scatter (%.4f): %.1f" % (SIGMA_EMP[0], abs(aen.mean() - 1.0) / SIGMA_EMP[0])) print() print("4. Ferrer-i-Cancho and Sole (2001) report two regimes, with the second") print(" regime steeper, near a = 2, beyond a crossover around rank 10^3 to 10^4.") w0 = np.array([r["win"][0] for r in rows if r["win"][0]]) w2 = np.array([r["win"][2] for r in rows if r["win"][2]]) print(" Club window r 1-100 : mean a = %.3f" % w0.mean()) print(" Club window r 1001-10000 : mean a = %.3f" % w2.mean()) print(" Same direction, same rough size.") return md, float(ben.mean()), float(aen.mean()) # --------------------------------------------------------------------------- # 12. Numbers the figures are drawn from. # --------------------------------------------------------------------------- def local_slopes(counts, nbins=16): V = counts.size edges = np.unique(np.round(np.logspace(0, math.log10(V), nbins + 1)).astype(np.int64)) xs, ys = [], [] for i in range(len(edges) - 1): lo, hi = int(edges[i]), int(edges[i + 1]) if hi <= lo: continue seg = counts[lo - 1:hi - 1] if hi > lo else counts[lo - 1:lo] if seg.size == 0: continue r = np.arange(lo, lo + seg.size, dtype=float) xs.append(float(np.exp(np.mean(np.log(r))))) ys.append(float(np.mean(seg))) xs, ys = np.array(xs), np.array(ys) rmid = np.sqrt(xs[1:] * xs[:-1]) slope = -np.diff(np.log(ys)) / np.diff(np.log(xs)) return rmid, slope def figure_data(rows, conv): head("FIGURE DATA. Every number the article's figures are drawn from") by = {r["label"]: r for r in rows} picks = ["Moby-Dick", "Pride and Prejudice", "Les miserables I", "Don Quijote", "Seitseman veljesta", "Divina Commedia"] print() print("FIG1 rank-frequency curves, rank then count, log-spaced ranks") for name in picks: r = by[name] c = r["counts"] ranks = np.unique(np.round(np.logspace(0, math.log10(c.size), 40)).astype(np.int64)) ranks = ranks[ranks <= c.size] print("FIG1 %s|%s|N=%d|V=%d" % (name, r["lang"], r["N"], r["V"])) print("FIG1pts " + " ".join("%d,%d" % (int(k), int(c[k - 1])) for k in ranks)) print() print("FIG2 per-book exponents: label|lang|year|N|V|a_zipf|se|a_zm|q|b|b_se|p_gof|hapax") for r in rows: print("FIG2 %s|%s|%d|%d|%d|%.4f|%.5f|%.4f|%.3f|%.4f|%.4f|%.3f|%.4f" % (r["label"], r["lang"], r["year"], r["N"], r["V"], r["a_z"], r["se_z"], r["a_m"], r["q"], r["b"], r["b_se"], r["p_gof"], r["hapax"])) print() print("FIG3 local slope against rank (binned), for three books") for name in ("Moby-Dick", "Don Quijote", "Seitseman veljesta"): rm, sl = local_slopes(by[name]["counts"]) print("FIG3 %s " % name + " ".join("%.1f,%.3f" % (a, b) for a, b in zip(rm, sl))) print() print("FIG4 Moby-Dick observed against both fitted forms, at log-spaced ranks") r = by["Moby-Dick"] c = r["counts"] V, N = r["V"], r["N"] rr = np.arange(1, V + 1, dtype=float) pz = rr ** (-r["a_z"]) pz = N * pz / pz.sum() pm = (rr + r["q"]) ** (-r["a_m"]) pm = N * pm / pm.sum() ranks = np.unique(np.round(np.logspace(0, math.log10(V), 44)).astype(np.int64)) ranks = ranks[ranks <= V] print("FIG4 params a_zipf=%.4f a_zm=%.4f q=%.3f N=%d V=%d" % (r["a_z"], r["a_m"], r["q"], N, V)) print("FIG4obs " + " ".join("%d,%d" % (k, int(c[k - 1])) for k in ranks)) print("FIG4zipf " + " ".join("%d,%.2f" % (k, pz[k - 1]) for k in ranks)) print("FIG4zm " + " ".join("%d,%.2f" % (k, pm[k - 1]) for k in ranks)) print() print("FIG5 validation convergence, from VALIDATION C at b_true = 1.95") for n in sorted(conv): m, sd, mse, bias = conv[n] print("FIG5 %d %.5f %.5f %.5f %+.5f" % (n, m, sd, mse, bias)) print() print("FIG5b tokeniser sensitivity and half-to-half scatter") print("FIG5b sigma_emp %.5f" % SIGMA_EMP[0]) print("FIG5b halfdiff " + " ".join("%.4f" % r["half_diff"] for r in rows)) # --------------------------------------------------------------------------- # 13. Main. # --------------------------------------------------------------------------- def main(): print("=" * 78) print("COUNTING EVERY WORD IN FIFTY BOOKS TO TEST A LAW NOBODY CAN EXPLAIN") print("Science Journaling Club, Volume 2 Issue 4, Summer 2026") print("=" * 78) print("run date : %s" % date.today().isoformat()) print("python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) print("master seed : %d" % MASTER_SEED) print("data : Project Gutenberg, https://www.gutenberg.org/") print("cache : analysis/data/gutenberg/") print() print("The club has no laboratory. This program is the experiment: it counts") print("words in fifty public-domain books and fits two families of distribution") print("to the counts by maximum likelihood. Nothing here was observed anywhere") print("but in those fifty files.") ss = np.random.SeedSequence(MASTER_SEED) sa, sb, sc, sd, se, sf = ss.spawn(6) validation_zeta() validation_sampler(sb) _, conv = validation_recovery(sc) validation_gof(sd) validation_rank(se) rows = load_corpus() t = time.time() for r in rows: r["tokens"] = tokenise(r["body"]) r["counts"], _ = counts_from_tokens(r["tokens"]) print() print("Tokenised %d books, %s tokens in total, in %.1f s." % (len(rows), format(sum(len(r["tokens"]) for r in rows), ","), time.time() - t)) quality_check(rows) analyse_books(rows, sf) head("RESULTS TABLE. Fifty books, one row each") print("%-32s %-5s %-6s %-9s %-7s %-7s %-7s %-6s %-6s %-6s %-5s %-6s %s" % ("title", "lang", "year", "tokens", "types", "a_zipf", "SE", "a_ZM", "q", "b", "xmin", "p_GOF", "hapax")) for r in rows: print("%-32s %-5s %-6d %-9s %-7s %-7.4f %-7.4f %-7.4f %-6.2f %-6.3f %-5d %-6.3f %.3f" % (r["label"][:32], r["lang"], r["year"], format(r["N"], ","), format(r["V"], ","), r["a_z"], r["se_z"], r["a_m"], r["q"], r["b"], r["xmin"], r["p_gof"], r["hapax"])) grp = group_tests(rows) tail_analysis(rows) mand = mandelbrot_analysis(rows) gofs = gof_summary(rows) hapax_check(rows) sigma_emp, sd_obs, Q2, df2 = split_half(rows) SIGMA_EMP[0] = sigma_emp published_comparison(rows) tokeniser_sensitivity(rows) figure_data(rows, conv) head("SUMMARY. The headline numbers") az = np.array([r["a_z"] for r in rows]) aen = np.array([r["a_z"] for r in rows if r["lang"] == "en"]) b = np.array([r["b"] for r in rows]) p = np.array([r["p_gof"] for r in rows]) print("books : %d" % len(rows)) print("written languages : %d" % len({r["lang"] for r in rows})) print("tokens counted : %s" % format(int(sum(r["N"] for r in rows)), ",")) print("word types counted : %s" % format(int(sum(r["V"] for r in rows)), ",")) print("rank exponent a, all books : %.4f, SD %.4f" % (az.mean(), az.std(ddof=1))) print("rank exponent a, English books : %.4f, SD %.4f" % (aen.mean(), aen.std(ddof=1))) print("range of a across the corpus : %.4f to %.4f" % (az.min(), az.max())) print("measured one-book repeatability : %.4f" % sigma_emp) print("spread of a in units of that : %.0f" % (az.std(ddof=1) / sigma_emp)) print("frequency exponent b, all books : %.3f, SD %.3f" % (b.mean(), b.std(ddof=1))) print("books whose pure power law survives GOF : %d of %d" % (int((p > 0.10).sum()), p.size)) print("books where Zipf-Mandelbrot wins on AIC : %d of %d" % (mand[0], len(rows))) print("Moby-Dick b, club : %.3f +/- %.3f" % ([r for r in rows if r["label"] == "Moby-Dick"][0]["b"], [r for r in rows if r["label"] == "Moby-Dick"][0]["b_se"])) print("Moby-Dick b, Clauset et al. 2009 : 1.95 +/- 0.02") print() print("total runtime : %.1f s" % (time.time() - T0)) print() rule() print("END OF OUTPUT") rule() if __name__ == "__main__": main()