Fitting both laws to real rank counts, the coverage curve, and the length trap that makes type token ratio misleading.
The notebook carries the source, the install step and every run below, so nothing has to be on your machine. Open it, then choose Runtime › Run all.
Colab opens it read only. Click Copy to Drive to keep
your changes. You can also download
empirical_laws.ipynb, or run it
locally: python3 empirical_laws.py.
python3 empirical_laws.py
The fits, the coverage table, and what happens to TTR when the text gets longer.
Download empirical_laws.py
· served verbatim at https://nlp.jcrlabz.com/code/empirical_laws.py
"""Zipf, Mandelbrot, Heaps and TTR, worked as Chapter 3 works them.
Four ideas, all of them a straight line hiding inside a curve.
facts the Brown totals every other section leans on
zipf fit alpha by least squares, and watch it move with the range
mandelbrot add the rank offset beta and watch the error drop
coverage how many tokens the top k words account for
heaps fit the vocabulary growth exponent, then extrapolate and fail
prep which cleaning step moves tokens, and which moves types
ttr why two type-token ratios are not comparable
python3 empirical_laws.py zipf
python3 empirical_laws.py all
The data is measured, not invented, and it is measured the way the course
grades it. Every Brown figure below comes from `autograder/nlp_pipeline.py`,
which is the Assignment 1 contract: its TOKEN_RE, lowercased, natural logs,
and fit_heaps sampling the growth curve every 1000 tokens. Run A1 against
your own corpus and the procedure will be the one shown here.
The budget-speech ranks are the fifteen printed on the Zipf slide.
Install: nothing, the Python standard library is enough
"""
import argparse
import math
# ------------------------------------------------------------------ data
# Indian budget speeches. (rank, frequency) as printed on the slide.
BUDGET = [(1, 59042), (2, 46110), (3, 35873), (4, 24661), (5, 22916),
(6, 14426), (7, 14366), (8, 11155), (9, 9760), (10, 9070),
(233, 479), (234, 478), (235, 477), (236, 476), (237, 476)]
BUDGET_V = 92366 # last rank in the frequency table
# Brown corpus, measured through autograder/nlp_pipeline.py.
BROWN_N = 1_023_734
BROWN_V = 44_324
BROWN_TOP = [("the", 70003), ("of", 36473), ("and", 28935), ("to", 26247),
("a", 23507), ("in", 21419), ("that", 10596), ("is", 10109),
("was", 9815), ("he", 9548)]
BROWN_HAPAX_PCT = 39.5 # types occurring exactly once
BROWN_LE2_PCT = 54.3 # types occurring twice or fewer
# alpha from fit_zipf over a rank range. The range is the point.
BROWN_ALPHA = [(100, 0.9816), (1000, 0.9668), (10000, 1.0726),
(BROWN_V, 1.3492)]
# (rank range, alpha, beta) from fit_mandelbrot
BROWN_MANDELBROT = [(1000, 0.9668, 0.00), (10000, 1.0984, 10.25),
(BROWN_V, 1.5180, 403.25)]
# cumulative share of all tokens held by the top k types, measured
BROWN_COVERAGE = [(10, 24.09), (100, 47.18), (1000, 69.13), (10000, 92.60)]
# vocabulary after N tokens, sampled off the same heaps_curve fit_heaps uses
BROWN_GROWTH = [(1_000, 468), (10_000, 2_584), (100_000, 12_962),
(500_000, 31_854), (1_000_000, 43_648)]
BROWN_HEAPS_K = 16.18
BROWN_HEAPS_B = 0.5764
# (genre, tokens, raw TTR, MSTTR over non-overlapping 8000-token windows)
BROWN_GENRES = [("press reportage", 90_089, 0.1359, 0.3077),
("humour", 18_450, 0.2552, 0.3021),
("general fiction", 58_605, 0.1458, 0.2744),
("romance", 58_966, 0.1307, 0.2520),
("learned", 164_036, 0.0885, 0.2478),
("government", 63_214, 0.1106, 0.2320)]
# effect of each preprocessing knob, as a share of the raw figures
BROWN_PREP = [("stop words removed", 55.5, 99.8),
("numbers stripped", 99.0, 97.5),
("Porter stemming", 100.0, 65.4),
("drop count <= 5", 94.1, 28.0)]
# ------------------------------------------------------------------ maths
def least_squares(xs, ys):
"""Slope and intercept of the line that minimises squared error.
Every law in this chapter is a power law, and every power law is a
straight line once you take logarithms of both axes. So the whole of
curve fitting here is this one function applied to logged data.
"""
n = len(xs)
mx, my = sum(xs) / n, sum(ys) / n
sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
sxx = sum((x - mx) ** 2 for x in xs)
slope = sxy / sxx
return slope, my - slope * mx
def sse(xs, ys, slope, intercept):
"""Sum of squared errors, in log space, of a fitted line."""
return sum((y - (slope * x + intercept)) ** 2 for x, y in zip(xs, ys))
def fit_zipf(points):
"""Fit f(r) = C / r^alpha. Returns (alpha, C, error)."""
xs = [math.log10(r) for r, _ in points]
ys = [math.log10(f) for _, f in points]
slope, intercept = least_squares(xs, ys)
return -slope, 10 ** intercept, sse(xs, ys, slope, intercept)
def fit_mandelbrot(points, beta_max=200.0, step=0.1):
"""Fit f(r) = C / (r + beta)^alpha by scanning beta.
beta is inside the logarithm, so it cannot be fitted by least squares.
We scan it and least-squares the rest. Coarse, and enough to show that
the valley is shallow.
"""
ys = [math.log10(f) for _, f in points]
best = None
beta = 0.0
while beta <= beta_max:
xs = [math.log10(r + beta) for r, _ in points]
slope, intercept = least_squares(xs, ys)
err = sse(xs, ys, slope, intercept)
if best is None or err < best[0]:
best = (err, -slope, beta, 10 ** intercept)
beta += step
err, alpha, beta, c = best
return alpha, beta, c, err
def harmonic(n, alpha):
"""H(n, alpha) = sum of r^-alpha for r = 1..n.
This is the normaliser that turns Zipf's law from a proportionality
into a probability distribution, and it is what makes coverage
computable.
"""
return sum(r ** (-alpha) for r in range(1, n + 1))
# ------------------------------------------------------------------ demos
def demo_facts():
print("\nBROWN the corpus the rest of this file measures\n")
print(f" tokens N {BROWN_N:>12,}")
print(f" types V {BROWN_V:>12,}")
print(f" TTR = V/N {BROWN_V/BROWN_N:>12.4f}")
top10 = sum(c for _, c in BROWN_TOP)
print(f"\n {'rank':>5}{'word':>10}{'count':>9}{'% of tokens':>13}")
print(" " + "-" * 37)
for i, (w, c) in enumerate(BROWN_TOP, 1):
print(f" {i:>5}{w:>10}{c:>9,}{100*c/BROWN_N:>12.2f}%")
print(f"\n The top ten types are {100*top10/BROWN_N:.2f}% of all tokens.")
print(f" {BROWN_HAPAX_PCT}% of the types occur exactly once, and "
f"{BROWN_LE2_PCT}% twice or fewer.")
print("\n These come from the A1 tokenizer. A stricter rule that keeps")
print(" only fully alphabetic tokens gives a different V and a different")
print(" hapax share. Neither is wrong. State the tokenizer, always.")
def demo_zipf():
print("\nZIPF f(r) = C / r^alpha, budget speeches\n")
# A student's first instinct: two points and a ratio.
f1 = dict(BUDGET)[1]
f10 = dict(BUDGET)[10]
alpha2 = math.log10(f1 / f10)
print(f" two-point estimate from r=1 and r=10")
print(f" alpha = log10({f1}/{f10}) = {alpha2:.4f}")
pred = f1 / 235 ** alpha2
print(f" predict f(235) = {f1} / 235^{alpha2:.4f} = {pred:.0f}"
f" actual 477 off by {100*(pred-477)/477:+.0f}%")
# The fix: use every point you have, in log space.
alpha, c, err = fit_zipf(BUDGET)
print(f"\n least squares over all {len(BUDGET)} points")
print(f" alpha = {alpha:.4f} C = {c:,.0f} SSE = {err:.4f}")
pred = c / 235 ** alpha
print(f" predict f(235) = {c:,.0f} / 235^{alpha:.4f} = {pred:.0f}"
f" actual 477 off by {100*(pred-477)/477:+.0f}%")
print("\n Two points threw away thirteen measurements and paid for it.")
print(" The head of the distribution is flatter than the body, so any")
print(" estimate read off the head alone understates alpha.")
# alpha is not a constant of the language. It is a constant of a range.
print(f"\n alpha depends on the rank range you fit, Brown corpus:")
for hi, a in BROWN_ALPHA:
print(f" ranks 1..{hi:<6} alpha = {a:.4f}")
print(" Report the range with the exponent, or the exponent means")
print(" nothing. The tail falls faster than the head, so the wider the")
print(" window, the steeper the fitted line.")
def demo_mandelbrot():
print("\nMANDELBROT f(r) = C / (r + beta)^alpha, budget speeches\n")
a0, c0, e0 = fit_zipf(BUDGET)
alpha, beta, c, err = fit_mandelbrot(BUDGET)
print(f" {'model':<14}{'alpha':>8}{'beta':>7}{'C':>12}{'SSE':>9}")
print(" " + "-" * 50)
print(f" {'Zipf':<14}{a0:>8.4f}{'0':>7}{c0:>12,.0f}{e0:>9.4f}")
print(f" {'Mandelbrot':<14}{alpha:>8.4f}{beta:>7.1f}{c:>12,.0f}{err:>9.4f}")
pred = c / (235 + beta) ** alpha
print(f"\n predict f(235) = {pred:.0f}, actual 477, "
f"off by {100*(pred-477)/477:+.0f}%")
print(f" The error drops by {100*(e0-err)/e0:.0f}% for one extra"
f" parameter.")
print("\n beta shifts every rank to the right before the power is")
print(" taken, which flattens the curve at small r and leaves the tail")
print(" alone. It is a head correction and nothing else.")
print("\n On Brown, beta depends entirely on how much of the tail you")
print(" fit, because beta is what absorbs the bend:")
print(f" {'ranks':<16}{'alpha':>8}{'beta':>10}")
for hi, a, b in BROWN_MANDELBROT:
print(f" 1..{hi:<13,}{a:>8.4f}{b:>10.2f}")
print(" Fit the head alone and beta is zero, so Mandelbrot buys nothing.")
print(" Fit every rank, which is what the A1 grader does, and beta runs")
print(" into the hundreds. A large beta is not a red flag. It is what")
print(" a single power law needs in order to cover a curve that bends.")
print("\n The valley is shallow. Nearby (alpha, beta) pairs fit almost")
print(" as well, so the pair is meaningful and neither number is:")
ys = [math.log10(f) for _, f in BUDGET]
for b in (0.0, 0.5, 1.0, 2.0, 5.0):
xs = [math.log10(r + b) for r, _ in BUDGET]
s, i = least_squares(xs, ys)
print(f" beta = {b:<5} best alpha = {-s:.4f} "
f"SSE = {sse(xs, ys, s, i):.4f}")
def demo_coverage():
print("\nCOVERAGE what fraction of tokens the top k types hold\n")
print(" Zipf as a probability, not a proportionality:")
print(" p(r) = r^-alpha / H(V, alpha), H(V,a) = sum r^-a, r=1..V")
print(" so the top k types cover H(k, alpha) / H(V, alpha).\n")
alpha = 0.9668 # Brown, fitted on ranks 1..1000
hv = harmonic(BROWN_V, alpha)
print(f" Brown, V = {BROWN_V:,}, alpha = {alpha} fitted on ranks 1..1000")
print(f" {'k':>7}{'model':>10}{'measured':>11}{'% of V':>9}")
print(" " + "-" * 37)
for k, measured in BROWN_COVERAGE:
model = 100 * harmonic(k, alpha) / hv
print(f" {k:>7,}{model:>9.2f}%{measured:>10.2f}%"
f"{100*k/BROWN_V:>8.2f}%")
print("\n Read the last column against the third. About 2% of the")
print(" vocabulary carries about 69% of the running text.")
print(" The model tracks the measurement at the head and drifts at the")
print(" tail, because a single alpha cannot bend twice.")
print("\n This is the number that justifies three later decisions.")
print(" a stop list is short because the head is short")
print(" a truncated vocabulary loses few tokens, not few types")
print(" a subword vocabulary of 32k spends its budget on the head")
def fit_heaps():
"""The (k, b) that autograder/nlp_pipeline.fit_heaps returns for Brown.
Not a re-fit of the five rows printed below. fit_heaps samples the
growth curve every 1000 tokens and fits all 1024 points, so quoting a
fit of five hand-picked marks would disagree with the grader.
"""
return BROWN_HEAPS_B, BROWN_HEAPS_K
def demo_heaps():
print("\nHEAPS V(N) = k N^b, Brown, sampled as fit_heaps samples it\n")
print(f" {'N':>11}{'V':>9}{'TTR':>9}")
print(" " + "-" * 29)
for n, v in BROWN_GROWTH:
print(f" {n:>11,}{v:>9,}{v/n:>9.4f}")
slope, k_all = fit_heaps()
print(f"\n fit_heaps over the whole 1024-point curve:")
print(f" b = {slope:.4f} k = {k_all:.2f}")
print(" b < 1, so vocabulary growth decelerates. It never stops,")
print(" because a power law has no asymptote.")
print("\n k and b are COUPLED, and this is the part most readers miss.")
print(" A fit that raises b must lower k to pass through the same data,")
print(" so a large b always arrives with a small k. Pin b on Brown and")
print(" the k that best fits moves with it:")
print(f" {'b pinned':>10}{'best k':>12}")
for pb, pk in [(0.4000, 155.93), (0.5000, 43.17), (0.6000, 11.95)]:
print(f" {pb:>10.4f}{pk:>12.2f}")
print(f" {slope:>10.4f}{k_all:>12.2f} (both fitted freely)")
print(" So quote the pair. 'k between 10 and 100 and b between 0.4 and")
print(" 0.6' is one statement about a curve, not two independent ones.")
# Extrapolation, which is where the law bites back.
early = [p for p in BROWN_GROWTH if p[0] <= 100_000]
xs_e = [math.log(n) for n, _ in early]
ys_e = [math.log(v) for _, v in early]
s_e, i_e = least_squares(xs_e, ys_e)
k_e = math.exp(i_e)
pred = k_e * 1_000_000 ** s_e
actual = dict(BROWN_GROWTH)[1_000_000]
print(f"\n Now fit only the first 100,000 tokens and extrapolate.")
print(f" On the full curve that prefix gives b = 0.7176, k = 3.50, and")
print(f" V(1,000,000) = 70,799 actual {actual:,} off by +62%")
print(" One decade of extrapolation, well over half again too many types.")
print(" b is not a constant of the language either. It drifts down as N")
print(" grows, because the log-log curve is gently concave.")
print(" Fit over the range you will actually operate in.")
# The link back to Zipf.
print(" Zipf and Heaps are one fact seen twice. If frequency falls as")
print(" r^-alpha with alpha > 1, then vocabulary grows as N^(1/alpha).")
alpha_all = BROWN_ALPHA[-1][1]
print(f" alpha over the full rank range = {alpha_all:.4f}")
print(f" so the predicted b = 1/alpha = {1/alpha_all:.4f}")
print(f" measured b = {slope:.4f}")
print(" Right order, wrong decimals, and the gap is honest. The")
print(" derivation assumes one alpha covers every rank, and the zipf")
print(" demo showed alpha moving from 0.98 to 1.35 across ranges.")
print(" Take the relationship, not the number: steeper frequency decay")
print(" means fewer rare words, which means slower vocabulary growth.")
print()
def demo_prep():
print("\nPREPROCESSING which knob moves tokens, which moves types\n")
print(f" {'operation':<22}{'tokens left':>13}{'types left':>12}")
print(" " + "-" * 47)
for name, tok_pct, typ_pct in BROWN_PREP:
print(f" {name:<22}{tok_pct:>12.1f}%{typ_pct:>11.1f}%")
print("\n Read the two columns against each other.")
print(" Stop-word removal takes 44% of the tokens and 0.2% of the")
print(" types. Zipf's law says exactly this: stop words are very few")
print(" types carrying very many tokens.")
print(" Porter stemming is the mirror image. It leaves every token in")
print(" place and collapses 35% of the vocabulary.")
print(" Dropping every type seen five times or fewer removes 72% of the")
print(" vocabulary and costs 5.9% of the tokens. That is the trade the")
print(" coverage table priced, and it is why a fixed subword vocabulary")
print(" of 32k is not the sacrifice it looks like.")
print("\n Preprocessing barely moves the Heaps fit: k goes 16.18 -> 17.33")
print(" and b goes 0.5764 -> 0.5958 once stop words and numbers go.")
print(" The curve is a property of the language, not of your cleaning.")
def demo_ttr():
print("\nTTR V/N, and why you may not compare two of them\n")
print(" Heaps gives V = K N^b, so TTR = V/N = K N^(b-1).")
print(" b < 1 makes the exponent negative, so TTR falls with length.")
print(" It is not a property of a text. It is a property of a text")
print(" at a length.\n")
b, k = fit_heaps() # the same fit the heaps demo prints
print(f" predicted from Brown's fit, k = {k:.2f}, b = {b:.4f}")
print(f" {'N':>11}{'TTR':>9}")
print(" " + "-" * 20)
for n in (1_000, 10_000, 100_000, 1_000_000):
print(f" {n:>11,}{k * n ** (b - 1):>9.4f}")
print(" Same author, same style, four different answers.\n")
print(" Six Brown genres. Raw TTR uses the whole genre. Standardised")
print(" TTR averages over non-overlapping 8,000-token windows.\n")
print(f" {'genre':<18}{'N':>9}{'raw TTR':>10}{'rank':>6}"
f"{'std TTR':>10}{'rank':>6}")
print(" " + "-" * 59)
raw_rank = {g[0]: i for i, g in
enumerate(sorted(BROWN_GENRES, key=lambda g: -g[2]), 1)}
std_rank = {g[0]: i for i, g in
enumerate(sorted(BROWN_GENRES, key=lambda g: -g[3]), 1)}
for name, n, raw, std in sorted(BROWN_GENRES, key=lambda g: -g[3]):
print(f" {name:<18}{n:>9,}{raw:>10.4f}{raw_rank[name]:>6}"
f"{std:>10.4f}{std_rank[name]:>6}")
print("\n Two rankings, two different answers.")
print(" Raw TTR calls humour the richest genre in Brown. Humour is")
print(" also the shortest genre in Brown, at 18,450 tokens against")
print(" 164,036 for learned prose. The measurement is reading length.")
print(" On equal windows press reportage comes first and humour")
print(" second, and press overtakes fiction, which raw TTR had the")
print(" other way round.")
print("\n Never compare type-token ratios across unequal samples.")
def main():
ap = argparse.ArgumentParser(description=__doc__)
ap.add_argument("which", nargs="?", default="all",
choices=["facts", "zipf", "mandelbrot", "coverage",
"heaps", "prep", "ttr", "all"])
a = ap.parse_args()
if a.which in ("facts", "all"):
demo_facts()
if a.which in ("zipf", "all"):
demo_zipf()
if a.which in ("mandelbrot", "all"):
demo_mandelbrot()
if a.which in ("coverage", "all"):
demo_coverage()
if a.which in ("heaps", "all"):
demo_heaps()
if a.which in ("prep", "all"):
demo_prep()
if a.which in ("ttr", "all"):
demo_ttr()
print()
if __name__ == "__main__":
main()