The chain behind Monte Carlo, PageRank and the next word
One 1906 idea, run four ways on this machine: the letters of Eugene Onegin, the neutrons at Los Alamos, the whole web, and a model writing its own text back to itself.
Probability
Markov Chains
Monte Carlo
PageRank
Python
Author
Ricky Macharm
Published
September 26, 2026
Three things you have almost certainly used this week run on one idea. One decided whether a wartime reactor would stay lit or run away. One decided which of the billions of pages on the web you were shown first. One decided the next word your phone offered you. The idea is a chain, and it is old enough to have been used on a poem.
A chain is a way of describing something that moves between states. A coin is heads or tails. A page is one of the pages on the web. A letter is a vowel or a consonant. What makes it a chain rather than a list of states is one restriction: where it goes next depends only on where it is now, and not on how it got there.
That sounds like a small restriction. It is the whole trick.
Put the states in a row, write down the chance of moving from each one to each other, and the future becomes arithmetic. You can draw a sample from it, average over it, or multiply a vector by it until the numbers stop moving. There is no need to remember history, because there is no history to remember.
The idea went into the literature through an argument about free will, and it reached the world’s computers through a card game played during a slow recovery from illness. Between those two points sits a counting exercise on a Russian poem that a mathematician did to win a fight. This post follows the whole line, and every number in it was produced by the code you will see, on this machine, a few minutes ago.
A fight about whether independence was necessary
The story starts with a theorem that had been standing for nearly two centuries. Jacob Bernoulli’s Ars Conjectandi was published in 1713, eight years after his death, and it contained the first serious proof of what we now call the law of large numbers: as you take more observations, the average of what you see closes in on the underlying probability.
Bernoulli’s proof assumed the observations were independent of each other. That assumption sat quietly under the result for two hundred years, the way a foundation sits under a house.
In the Russia of the early 1900s, one mathematician decided to make the assumption do political work. Pavel Nekrasov argued from statistical regularities in social data to a conclusion about human freedom: if births, marriages and crimes settle into steady averages, he reasoned, then individual human acts must be independent of one another, and that independence pointed to something beyond material cause.
One of the tables behind the argument was Adolphe Quetelet’s Belgian returns. Over five consecutive years, 1841 to 1845, the number of marriages registered in Belgium was 29,876, 29,023, 28,220, 29,326 and 29,210. The whole span stays inside 5.7% of its own average, and the typical year moves 2.81%. Nobody arranged that. A country of a few million people simply kept arriving at about twenty-nine thousand marriages a year.
Quetelet’s own image for regularity was not about weddings. Reading an 1838 essay on French moral statistics, he stopped at a line about crime: “there is a budget which we pay with frightful regularity — it is that of prisons, dungeons, and scaffolds.” He read it with surprise, then took it up a page later, calling the budget of the scaffold and the prisons more regularly paid than the financial one.
The same counting is still running. Germany recorded 387,423 marriages in 2012 and 390,743 in 2022 — eleven years apart, and 0.86% apart — while its marriage rate stayed between 4.3 and 5.4 per 1,000 residents across the ten years to 2022, averaging 4.8. Births in the United States never left the band 3,641,990 to 4,232,701 over the twenty-four years to 2023, and the typical year moved 1.15% from the one before.
"""The same regularity, counted in Belgium in the 1840s and in two countries today."""import csvimport statistics as strows = [r for r in csv.reader(open("national_regularity.csv", encoding="utf-8")) if r andnot r[0].startswith("#")]series = {}for name, country, year, value in rows[1:]: series.setdefault(name, []).append((int(year), float(value)))def band(name, fmt="{:,.0f}"): pts =sorted(series[name]) years = [y for y, _ in pts] vals = [v for _, v in pts] moves = [100*abs(vals[i +1] - vals[i]) / vals[i] for i inrange(len(vals) -1)]print(f"{name}{min(years)}-{max(years)}, {len(vals)} years")print(f" lowest {fmt.format(min(vals))} highest {fmt.format(max(vals))} average {fmt.format(st.mean(vals))}")print(f" typical move year to year {st.median(moves):.2f}% biggest {max(moves):.2f}%")print(f" the whole span sits inside {100* (max(vals) -min(vals)) / st.mean(vals):.1f}% of the average")return valsprint("the numbers Nekrasov was arguing about are still being counted.")print()print("Belgian marriages, the five years Quetelet tabulated")be = band("be_marriages")print(" "+" ".join(f"{int(v):,}"for v in be))print()print("German marriages, a century and a half later")de = band("de_marriages")print(f" {int(de[0]):,} in 2012 against {int(de[-1]):,} in 2022 — {100*abs(de[-1] - de[0]) / de[0]:.2f}% apart, eleven years on")print()print("Births in the United States")us = band("us_births")print()rate = [v for _, v insorted(series["de_marriages_per_1000"])]print(f"and the German marriage rate never left {min(rate):.1f} to {max(rate):.1f} marriages per 1,000 residents")print(f"in the ten years to 2022, averaging {st.mean(rate):.1f}.")
the numbers Nekrasov was arguing about are still being counted.
Belgian marriages, the five years Quetelet tabulated
be_marriages 1841-1845, 5 years
lowest 28,220 highest 29,876 average 29,131
typical move year to year 2.81% biggest 3.92%
the whole span sits inside 5.7% of the average
29,876 29,023 28,220 29,326 29,210
German marriages, a century and a half later
de_marriages 2012-2022, 11 years
lowest 357,785 highest 449,466 average 395,696
typical move year to year 3.91% biggest 10.33%
the whole span sits inside 23.2% of the average
387,423 in 2012 against 390,743 in 2022 — 0.86% apart, eleven years on
Births in the United States
us_births 2000-2023, 24 years
lowest 3,641,990 highest 4,232,701 average 3,935,714
typical move year to year 1.15% biggest 3.46%
the whole span sits inside 15.0% of the average
and the German marriage rate never left 4.3 to 5.4 marriages per 1,000 residents
in the ten years to 2022, averaging 4.8.
The historian of probability Eugene Seneta has traced this argument in detail. What Nekrasov was doing, in effect, was reading a theorem backwards. The law of large numbers says that independence produces steady averages. Nekrasov wanted it to say that steady averages prove independence.
None of that proves the individual acts were independent. The yearly totals wander — one German year moved 10% — and a chain behaves exactly this way: because each year starts from the last, the aggregate moves slowly while the decisions inside it stay free. Nekrasov needed the steadiness to be evidence of independence. It is evidence of arithmetic.
Andrey Markov would not let it stand. He called this use of probability an abuse of mathematics, and rather than argue the philosophy he went after the implication itself: if steady averages require independence, then he would produce steady averages from something that was plainly not independent.
He chose the letters of a poem.
Twenty thousand letters of Eugene Onegin
Markov took the first chapter of Pushkin’s Eugene Onegin, sixteen stanzas of the second chapter, and stripped it to bare letters — no punctuation, no spaces, and none of the Russian hard and soft signs. That left him 20,000 characters to work with.
He sorted the letters into two kinds, vowels and consonants, and counted. The result, published in 1913 in the Bulletin of the Imperial Academy of Sciences of St Petersburg, was 43.19% vowels. If letters were independent of one another, then a vowel would be followed by a vowel about 0.4319 squared of the time, which is 18.65%. Markov counted 1,104 vowel-vowel pairs among his 20,000 letters, or 5.52%.
That gap is the argument. Letters are not independent — Russian vowels and consonants alternate far more than chance would allow — and yet the average still settled down. Independence was not required.
Here is that count, redone on a public-domain text of the poem. The code reads the letters, drops the hard and soft signs, and takes the first 20,000.
"""Markov's count, redone: the letters of Eugene Onegin."""text =open("onegin_markov_scope.txt", encoding="utf-8").read()# Markov counted letters only: no punctuation, no spaces, and he left out the# Russian hard and soft signs. Everything here is lowercase Cyrillic.def letters(s, drop="ъь"):return"".join(c for c in s.lower() if"а"<= c <="я"and c notin drop)VOWELS =set("аеёиоуыэюя") | {"й"} # as it turns out, he counted й as a vowelL = letters(text)[:20000]n_v =sum(1for c in L if c in VOWELS)vv =sum(1for a, b inzip(L, L[1:]) if a in VOWELS and b in VOWELS)p = n_v /len(L)print(f"letters counted {len(L):,}")print(f"vowels {n_v:,} ({p:.2%})")print(f"vowel after vowel {vv:,} ({vv / (len(L) -1):.2%})")print(f"if letters were independent {p **2:.2%} <- what independence predicts")print(f"p(vowel | vowel) {vv / n_v:.4f}")print("\nMarkov published 1913: 43.19% vowels | 5.52% vowel-vowel | p = 0.128")
letters counted 20,000
vowels 8,640 (43.20%)
vowel after vowel 1,110 (5.55%)
if letters were independent 18.66% <- what independence predicts
p(vowel | vowel) 0.1285
Markov published 1913: 43.19% vowels | 5.52% vowel-vowel | p = 0.128
That “as it turns out” in the code is doing real work, so it is worth stopping on. Markov’s paper does not spell out which letters he counted as vowels, and a different edition of the poem is not the same text he held in his hands. What the paper does give is three numbers — 43.19%, 5.52%, and the conditional probability of 0.128 — and those three are enough to work out how he was classifying letters, because only one classification reproduces all three.
"""Which letters did he call vowels? Only one convention reproduces his three figures."""print(f"{'variant':<44}{'vowels':>9}{'vowel-vowel':>13}{'p(v|v)':>9}")for label, drop, extra in [ ("everything kept", "", set()), ("hard and soft signs dropped (ъ, ь)", "ъь", set()), ("that, and й counted as a vowel", "ъь", {"й"}),]: L2 = letters(text, drop)[:20000] V =set("аеёиоуыэюя") | extra nv =sum(1for c in L2 if c in V) vv2 =sum(1for a, b inzip(L2, L2[1:]) if a in V and b in V)print(f"{label:<44}{nv /len(L2):>8.2%}{vv2 / (len(L2) -1):>12.2%}{vv2 / nv:>9.3f}")print(f"{'Markov, published 1913':<44}{0.4319:>8.2%}{0.0552:>12.2%}{0.128:>9.3f}")
variant vowels vowel-vowel p(v|v)
everything kept 40.28% 3.02% 0.075
hard and soft signs dropped (ъ, ь) 41.16% 3.11% 0.075
that, and й counted as a vowel 43.20% 5.55% 0.128
Markov, published 1913 43.19% 5.52% 0.128
Keep every letter and you get 40.28% vowels and vowel-vowel pairs at 3.02% — too low on both. Drop the hard and soft signs as he did and it rises to 41.16% and 3.11%, still short. Drop them and count й as a vowel and the three figures land at 43.20%, 5.55% and 0.1285 against his published 43.19%, 5.52% and 0.128.
Reproducing a dead man’s arithmetic within three thousandths, on a Soviet edition of the poem printed forty-six years after his paper, is how you find out what he meant.
Source: A. A. Markov, “An Example of Statistical Investigation of the Text Eugene Onegin Concerning the Connection of Samples in Chains,” 1913, English translation in Science in Context 19(4), 2006, 591–600; Eugene Seneta, “Statistical Regularity and Free Will: L. A. J. Quetelet and P. A. Nekrasov,” International Statistical Review 71, 2003, 319–334.
Figure 1: The two-state chain, built from Markov’s own counts: a vowel follows a vowel 12.85% of the time, a consonant follows a consonant 33.71%.
The chain, not the average
The move that survived is the one in the next block. Instead of counting the letters, Markov wrote down the chance of each letter following each state: after a vowel, a vowel came 12.85% of the time; after a consonant, a vowel came 66.29% of the time.
Then he generated letters from those two numbers alone. The poem was never consulted again, and the averages still settled down.
"""The same letters as a two-state chain, generating text from two numbers."""import numpy as npstates = ["v"if c in VOWELS else"c"for c in L]counts = {"v": {"v": 0, "c": 0}, "c": {"v": 0, "c": 0}}for a, b inzip(states, states[1:]): counts[a][b] +=1P = {s: {t: counts[s][t] /sum(counts[s].values()) for t in ("v", "c")} for s in ("v", "c")}print("transition probabilities from the poem itself")for s in ("v", "c"):print(f" after a {'vowel'if s =='v'else'consonant'}: "f"vowel {P[s]['v']:.4f} | consonant {P[s]['c']:.4f}")rng = np.random.default_rng(1913)start = states[0]shares = []for run inrange(20): state, out = start, []for _ inrange(20000): state ="v"if rng.random() < P[state]["v"] else"c" out.append(state) shares.append(out.count("v") /len(out))print(f"\n20 runs of 20,000 letters from the chain alone: {np.mean(shares):.2%} vowels on average")print(f"the poem itself: {states.count('v') /len(states):.2%} vowels")
transition probabilities from the poem itself
after a vowel: vowel 0.1285 | consonant 0.8715
after a consonant: vowel 0.6629 | consonant 0.3371
20 runs of 20,000 letters from the chain alone: 43.07% vowels on average
the poem itself: 43.20% vowels
Twenty runs of twenty thousand letters, drawn from two numbers, average 43.07% vowels. The poem itself is 43.20%. Nekrasov’s implication is dead: the regularity survives the loss of independence, and it survives it completely.
Markov’s 1906 paper had already carried the theorem where he needed it. Its title translates as Extension of the Law of Large Numbers to Quantities Depending on Each Other, and it is the ancestor of everything that follows in this post. What he did with the poem in 1913 was demonstrate the result on something nobody could accuse of being ordered.
What the restriction buys
A chain is memoryless, and that is the reason it is used. To predict the next letter you do not need the poem so far, only the kind of letter you are standing on. The state is the memory.
That property is what makes the arithmetic cheap. A system with a thousand states needs a thousand-by-thousand table of moves, and nothing else. You can multiply a vector by that table and watch the distribution settle; you can draw paths through it; you can count how often each state is visited in the long run. All of it is multiplication and sampling.
The cost is that the chain forgets on purpose. A chain cannot tell you that a long run of vowels is overdue for a consonant, or that a sentence has a subject waiting for a verb, because none of that is in the state. Everything that matters has to be squeezed into “where am I right now”.
Which is exactly why the restriction is useful, and exactly why it eventually breaks. Both halves of that sentence get their turn below. First the useful half, in a desert in New Mexico. ## A solitaire game in Los Alamos
In early 1946, Stanislaw Ulam was in a hospital bed in Los Alamos with encephalitis, and he was playing a great deal of solitaire. He was also meant to be thinking about how neutrons behave inside a piece of plutonium, which is a hard problem in a way that is worth stating plainly: a neutron’s life is a sequence of collisions, and each collision depends on the energy and direction the neutron arrived with.
Ulam’s observation was that you did not have to solve that sequence. You could play it. Deal a neutron a random collision, then another, and follow it until it is absorbed, leaks out, or splits an atom. Do that a few thousand times and count what happened.
John von Neumann, who had the computing machinery in mind, worked the idea up with him. The two of them noted that this is a chain: the neutron’s next collision depends on its current state and nothing earlier. Their 1949 paper in the Journal of the American Statistical Association describes cascading nuclear processes as Markoff chains and follows neutron transport as a sequence of random histories.
The name came from Nicholas Metropolis, and it came from gambling. Ulam had an uncle who gambled away his money at the casino in Monte Carlo, and the method was a game of chance played out at scale.
The number Ulam picked as his example of a game too large to enumerate was the number of ways to order a 52-card deck: 52 factorial, about 8.07 followed by 67 zeros. That figure is often given as the number of solitaire games, which is not quite it — it is the number of arrangements of the deck, and it is already a number no one will ever list.
The problem the method was built to answer has a name. Call it k: the average number of neutrons that the next generation of fissions produces, per neutron in this one. The Nuclear Regulatory Commission states the three cases in its reactor physics manual, and they are the whole of reactor design:
k below 1: the population shrinks and the reaction dies. k equals 1: the population holds steady and the reaction is self-sustaining. k above 1: the population grows, generation after generation.
The first Monte Carlo calculations ran on ENIAC, the Electronic Numerical Integrator and Computer, in the spring of 1948. By December of that year the same machine was being used on reactor problems for Argonne National Laboratory. The method that began as a way to follow one neutron became the standard way to design a core.
You can build the same thing in a few dozen lines. Each neutron travels, scatters, is absorbed, leaks out, or splits an atom; one knob sets how often it splits. Everything else is arithmetic. Here is k measured for seven settings of that knob, forty thousand neutrons each.
"""Neutrons as a chain: the multiplication factor k, measured."""import numpy as npSTATE = ["travel", "scatter", "absorbed", "left", "fission"]def table(fission_chance, fast=False):"""The five outcomes, scaled so the fission share is what we set.""" base = {"scatter": 0.40, "absorbed": 0.10, "left": 0.20, "fission": 0.30}if fast: base = {"scatter": 0.30, "absorbed": 0.20, "left": 0.30, "fission": 0.20} rest =1.0- fission_chance w = {k: v for k, v in base.items() if k !="fission"} scale = rest /sum(w.values()) row = np.zeros(len(STATE))for k, v in w.items(): row[STATE.index(k)] = v * scale row[STATE.index("fission")] = fission_chance t = np.tile(row, (len(STATE), 1))for dead in ("absorbed", "left"): # a neutron that is gone stays gone t[STATE.index(dead)] =0.0 t[STATE.index(dead), STATE.index(dead)] =1.0return tdef simulate(rng, probs, n_neutrons):"""Follow each neutron until it is absorbed, lost, or splits an atom.""" alive = np.ones(n_neutrons, dtype=bool) state = np.zeros(n_neutrons, dtype=np.int8) # all start in travel released = np.zeros(n_neutrons, dtype=np.int8) # what each neutron freesfor _ inrange(400): idx = np.flatnonzero(alive)if idx.size ==0:break u = rng.random(idx.size) nxt = (u[:, None] > np.cumsum(probs[state[idx]], axis=1)).sum(axis=1) state[idx] = nxt f = idx[nxt ==4] released[f] += rng.integers(2, 4, size=f.size) # a split frees 2 or 3 state[f] =0# and the neutron walks on alive[idx[(nxt ==2) | (nxt ==3)]] =False# absorbed or leaked outreturn releasedrng = np.random.default_rng(1946)print("one knob: the chance a neutron splits an atom instead of leaving")print(f"{'fission chance':>15}{'k (neutrons per neutron)':>27}")for chance in (0.05, 0.10, 0.145, 0.15, 0.20, 0.25, 0.30):print(f"{chance:>15.3f}{simulate(rng, table(chance), 40000).mean():>27.3f}")
one knob: the chance a neutron splits an atom instead of leaving
fission chance k (neutrons per neutron)
0.050 0.307
0.100 0.657
0.145 0.982
0.150 1.031
0.200 1.476
0.250 1.983
0.300 2.527
The knob moves k across the critical line: at a fission chance of 0.145 the core measures 0.982, at 0.15 it measures 1.031. Somewhere between those two settings sits the difference between a design that dies and a design that holds.
Figure 2: Reactivity against the one knob. The curve crosses k = 1 between a fission chance of 0.10 and 0.15.
Plot that against a population of neutrons and the three regimes stop being words and become a table. Four hundred neutrons, eight generations.
"""The same core, three settings, run generation by generation."""print("and what that does to a population of 400 neutrons, over 8 generations")CORES = (0.05, 0.145, 0.30)drawn = {c: simulate(np.random.default_rng(7), table(c), 40000) for c in CORES}ks = {c: float(drawn[c].mean()) for c in CORES}print(f"{'generation':>11}"+"".join(f"{'k = '+format(ks[c], '.2f'):>12}"for c in CORES))pop = {c: 400for c in CORES}print(f"{0:>11}"+"".join(f"{pop[c]:>12,}"for c in CORES))rngs = {c: np.random.default_rng(int(c *1000)) for c in CORES}for gen inrange(1, 9): row = []for c in CORES: pop[c] =int(np.sum(rngs[c].choice(drawn[c], size=pop[c]))) row.append(pop[c])print(f"{gen:>11}"+"".join(f"{x:>12,}"for x in row))
and what that does to a population of 400 neutrons, over 8 generations
generation k = 0.30 k = 0.99 k = 2.52
0 400 400 400
1 118 420 1,015
2 41 395 2,478
3 7 355 6,367
4 5 355 16,028
5 0 308 40,379
6 0 268 101,438
7 0 281 256,613
8 0 258 646,325
Four hundred neutrons at k = 0.30 are five and then nothing. At k = 2.52 they are 646,325 in eight generations. In between, at k = 0.99, they wander — up to 420, down to 258 — and that wandering is the point of the last measurement.
A core whose k sits at 1 does not behave like a core whose k is 1. Run it three hundred times and the average lands near 1 while individual runs do anything.
"""The same core, 300 times: why one run proves nothing."""rng = np.random.default_rng(99)runs = np.array([simulate(rng, table(0.20, fast=True), 300).mean() for _ inrange(300)])print(f"300 separate runs of the same core (fast neutrons):")print(f" mean k {runs.mean():.3f} | smallest {runs.min():.3f} | largest {runs.max():.3f}")print(f" runs that died out (k < 1): {(runs <1).mean():.0%} of them")
300 separate runs of the same core (fast neutrons):
mean k 0.997 | smallest 0.650 | largest 1.337
runs that died out (k < 1): 53% of them
Three hundred runs of one core, mean k = 0.997. More than half of them — 53% — finished below 1 and died. The Los Alamos group was not building one reactor; it was sampling a distribution of them, which is precisely the thing a chain simulation is for.
A note on the Trinity test, since that is where this arithmetic was first pointed. The device detonated on 16 July 1945 was a plutonium bomb, and the approximately 6 kilograms usually quoted for it refer to the plutonium content of the core, not the weight of the bomb. Published yields are lower than the roundest retelling: the Department of Energy’s history gives about 21 kilotons, and an Air Force history gives 18.6. The frequently repeated “nearly 25,000 tons of TNT” is on the high side of the official figures.
Source: Terrence R. Fehner and F. G. Gosling, Origins of the Nevada Test Site, U.S. Department of Energy, 2000; Beck et al., “Accounting for Unfissioned Plutonium from the Trinity Atomic Bomb Test,” Health Physics 119, 2020.
The web as a chain
By the mid-1990s the web had a search problem. Pages were being matched to queries mostly by the words on them, and the words on a page are chosen by whoever wrote the page, so a page could be pushed up the list by stuffing it with terms. Sergey Brin and Larry Page, both then doctoral candidates at Stanford, proposed in 1998 that a page’s importance should be decided by who links to it.
A link was a vote, in their formulation, and a page’s vote was divided evenly among the pages it linked to. Two pages linking to you are worth more than one, and a single link from a page that everyone else links to is worth more than a hundred links from nowhere.
The result is a number per page, and the way to compute it is a chain. Imagine a reader who clicks links forever at random, and ask what fraction of the time they spend on each page. That fraction is the page’s rank. Brin and Page wrote it as the principal eigenvector of the normalized link matrix and set the damping factor — the chance the reader follows a link rather than jumping to a random page — at 0.85.
Here is a toy web of four pages, scored by repeatedly multiplying a vector by the table of moves.
"""The web as a chain: four pages, scored by repeated multiplication."""import numpy as npdef rank(edges, pages, d=0.85, iters=100):"""edges: {page: [pages it links to]}. d: chance the surfer follows a link.""" n =len(pages) idx = {p: i for i, p inenumerate(pages)} M = np.zeros((n, n))for src, outs in edges.items():if outs:for dst in outs: M[idx[dst], idx[src]] =1/len(outs) # votes split evenly dangling = [idx[p] for p in pages ifnot edges.get(p)] v = np.ones(n) / nifnot d:for _ inrange(iters): v = M @ vreturn v, v.sum() # no jump: mass drains through the dead endfor d0 in dangling: # with the jump on, a dead end sends the surfer anywhere M[:, d0] =1/ nfor _ inrange(iters): v = d * (M @ v) + (1- d) / nreturn v, v.sum()pages = ["Amy", "Ben", "Chris", "Dan"]web = {"Amy": ["Ben"], "Ben": ["Amy", "Chris", "Dan"], "Chris": ["Amy", "Dan"], "Dan": ["Ben"]}v, total = rank(web, pages)print("a toy web: Amy -> Ben; Ben -> Amy, Chris, Dan; Chris -> Amy, Dan; Dan -> Ben")for p, s insorted(zip(pages, v), key=lambda t: -t[1]):print(f" {p:<6}{s:.4f}")print(f" total mass {total:.4f}")
a toy web: Amy -> Ben; Ben -> Amy, Chris, Dan; Chris -> Amy, Dan; Dan -> Ben
Ben 0.4092
Amy 0.2187
Dan 0.2187
Chris 0.1534
total mass 1.0000
Ben takes 0.41 of the reader’s time, Amy and Dan 0.22 each, Chris 0.15. Ben wins because the links into him come from pages that are themselves linked to, while Chris’s single vote has to be split two ways.
Figure 3: The toy web, ranked. Four pages, four scores, one random reader.
Now the failure that forced the damping factor into the design. Give Dan no links at all — a page with nothing to click — and run the same arithmetic without letting the reader jump anywhere.
"""A page with no links: what the random jump is actually for."""broken = {"Amy": ["Ben"], "Ben": ["Amy", "Chris", "Dan"], "Chris": ["Amy", "Dan"], "Dan": []}v0, t0 = rank(broken, pages, d=0)v1, t1 = rank(broken, pages)print(f" follow links only (d = 0): total mass {t0:.3f}")print(f" with the 15% random jump (d = 0.85): total mass {t1:.3f}")
follow links only (d = 0): total mass 0.000
with the 15% random jump (d = 0.85): total mass 1.000
The mass goes to zero. A reader who follows links and finds none has nowhere to go, so the walk dies and the ranking has nothing left to rank. The 15% jump is what keeps the walk alive: it teleports the reader to a random page, which is also why the rank of a page nobody links to is small but never zero.
Which raises the question every site owner asked in the years after 1998: if links are votes, why not build pages that vote for you? Hide a hundred pages that link only to your page, and see what it buys.
"""A hundred pages linking only to Amy."""spam = [f"spam{i:03d}"for i inrange(100)]big =dict(web)for s in spam: big[s] = ["Amy"]bv, _ = rank(big, pages + spam)print(f" Amy {v[0]:.4f} -> {bv[pages.index('Amy')]:.4f} (Ben still {bv[pages.index('Ben')]:.4f})")print(f" the best-placed spam page: {bv[len(pages):].max():.4f}, and all 100 together hold "f"{bv[len(pages):].sum():.3f} of the vote")
Amy 0.2187 -> 0.2652 (Ben still 0.3480)
the best-placed spam page: 0.0014, and all 100 together hold 0.144 of the vote
Amy rises from 0.2187 to 0.2652 and Ben stays on top at 0.3480. The hundred fake pages hold 0.144 of the total between them, and the best of them is worth 0.0014 — because each one has no links pointing at it, so under this arithmetic there is nothing to inherit.
That is the shape of the defence, and it is worth being honest about the shape of the attack. A link farm does move the number; it just cannot buy the top spot, and the pages doing the work stay at the bottom where nobody will find them.
Google incorporated in August 1998, after a rename from BackRub. The name is a misspelling of googol, the word for ten to the hundredth power, and Brin and Page’s own paper calls Google “a common spelling of googol”. The more specific story — that a fellow student typed the misspelling while registering the domain — comes from a Stanford account written years later, so treat it as a good story rather than a document.
Alphabet was worth about two trillion dollars in July 2025 and more than four trillion by early 2026. Yahoo, for the record, was founded in 1994 by Jerry Yang and David Filo, and Masayoshi Son did offer roughly 105 million dollars for a third of it the following year.
Source: Sergey Brin and Lawrence Page, “The Anatomy of a Large-Scale Hypertextual Web Search Engine,” Computer Networks and ISDN Systems 30, 1998, 107–117; Randall Lane, “Master of the Internet,” Forbes, 5 July 1999; CNBC, January 2026.
Shannon’s approximation, and what happens when a model reads itself
Claude Shannon’s 1948 paper A Mathematical Theory of Communication is where information theory starts, and it contains a demonstration that belongs in this post: he approximated English with chains and printed the results to show how much structure each order adds.
At the first order, letters are drawn independently and the output is noise. At higher orders, the chain conditions on the previous few letters or words, and the text starts to look like language. Here is the same experiment on a public-domain novel, with the information content measured at each order.
"""Shannon's approximations to English, with bits per token measured."""import reimport numpy as npfrom collections import Counter, defaultdictbook =open("corpus_pride_and_prejudice.txt", encoding="utf-8").read()words = re.findall(r"[a-z']+", book.lower())chars = re.findall(r"[a-z ]", book.lower())print(f"corpus: {len(words):,} words of Pride and Prejudice (Project Gutenberg ebook 1342)")def chain_of(tokens, order): c = defaultdict(Counter)for i inrange(len(tokens) - order): c[tuple(tokens[i:i + order])][tokens[i + order]] +=1return cdef bits(chain): total =sum(sum(v.values()) for v in chain.values()) h =0.0for counter in chain.values(): n =sum(counter.values()) p = np.array(list(counter.values())) / n h += n / total *float(-(p * np.log2(p)).sum())return hdef speak(chain, order, rng, n, tokens): i =int(rng.integers(0, len(tokens) - order -1)) state, out =tuple(tokens[i:i + order]), [] out.extend(state)for _ inrange(n):if state notin chain:break options, weights =zip(*chain[state].items()) state = (state + (rng.choice(options, p=np.array(weights) /sum(weights)),))[1:] out.append(state[-1])return outrng = np.random.default_rng(1948)print("\nword by word, learning from the book:")for order in (1, 2, 3, 4): c = chain_of(words, order)print(f"\n order {order} ({bits(c):5.2f} bits per word, {len(c):,} states)")print(f" \"{' '.join(speak(c, order, rng, 26, words))[:200]}\"")
corpus: 43,352 words of Pride and Prejudice (Project Gutenberg ebook 1342)
word by word, learning from the book:
order 1 ( 4.74 bits per word, 4,395 states)
"gave her release from the subject to contrast to judge very little unwillingness but darcy that my dear mr darcy was given her manners gave a sunday"
order 2 ( 1.36 bits per word, 24,870 states)
"son but this word cynical is one of humanity and especially to her daughter settled at netherfield said mrs bennet elizabeth and it is no sacrifice to join"
order 3 ( 0.21 bits per word, 39,044 states)
"in silence but was not convinced their behaviour at the assembly had not been gone long before it rained hard her sisters were uneasy for her but her mother"
order 4 ( 0.03 bits per word, 42,487 states)
"know to which of his fair cousins the excellence of its cookery was owing but here he was set right by mrs bennet who assured him with some asperity that"
The numbers fall as the order rises: 4.74 bits per word at order 1, 1.36 at order 2, 0.21 at order 3, 0.03 at order 4. That is Shannon’s point restated for words — the more context you allow the chain, the less surprise is left in each new word.
Do not read the low numbers as a good model. The measure collapses because the corpus is far too small to fill the higher orders, so most states are seen once and their continuation is a single recorded word. The chain is not learning the language, it is reciting.
That shows up when you ask a blunter question: of the word sequences this thing writes, how many were ever actually written, in that order, in the book?
"""How much of what the chain writes was ever actually written?"""for order in (2, 3, 4): c = chain_of(words, order) real = {tuple(words[i:i + order +1]) for i inrange(len(words) - order)} text = speak(c, order, rng, 300, words) hit =sum(1for i inrange(len(text) - order) iftuple(text[i:i + order +1]) in real) single =sum(1for v in c.values() iflen(v) ==1) /len(c)print(f" order {order}: {hit}/{len(text) - order} generated {order +1}-word runs appear "f"in the book ({hit / (len(text) - order):.1%}); {single:.1%} of its states have only "f"one recorded continuation")
order 2: 300/300 generated 3-word runs appear in the book (100.0%); 81.3% of its states have only one recorded continuation
order 3: 300/300 generated 4-word runs appear in the book (100.0%); 94.6% of its states have only one recorded continuation
order 4: 300/300 generated 5-word runs appear in the book (100.0%); 98.8% of its states have only one recorded continuation
Every one of them, at every order. With 43,000 words to learn from, 81% of the two-word states and 99% of the four-word states have exactly one continuation on record, so the chain mostly walks a path the book already contains. What looks like fluent text is a quotation engine with small steps.
Shannon’s own examples make the same point more elegantly. His third-order letter approximation, printed in the 1948 paper, reads: “IN NO IST LAT WHEY CRATICT FROURE BIRS GROCID PONDENOME OF DEMONSTURES OF THE REPTAGIN IS REGOACTIONA OF CRE”.
Modern language models are the same idea with two changes: the states are built from tokens rather than letters, and instead of conditioning on a fixed window of the last few, attention weighs every earlier token in the text. That last change is the whole difference between a chain and a transformer, and it is what lets a model use a noun from four paragraphs ago.
Now the part that has no precedent in 1913. Take the chain trained on the book, let it write, then train a new chain on what it wrote, then a new one on that, five times over. This is what happened when the web started filling up with model output.
"""Feed each chain its own output, five generations over, eight chains at once."""import statisticsprint("feed the chain its own output, five generations over, eight chains at once:")runs, sample = [], ""for seed inrange(8): r = np.random.default_rng(1000+ seed) text, row = words, []for gen inrange(1, 6): c = chain_of(text, 2) text = speak(c, 2, r, 4000, words if gen ==1else text) row.append((len(set(text)), len(c)))if seed ==0and gen ==5: sample =" ".join(text[100:150]) runs.append(row)print(f"{'generation':>11}{'vocabulary':>14}{'distinct pairs':>16}")for gen inrange(5): v = statistics.mean(r[gen][0] for r in runs) p = statistics.mean(r[gen][1] for r in runs)print(f"{gen +1:>11}{v:>14,.0f}{p:>16,.0f}")print(f"\nand what it writes by then: \"{sample}\"")
feed the chain its own output, five generations over, eight chains at once:
generation vocabulary distinct pairs
1 1,070 24,870
2 676 2,995
3 403 1,613
4 325 908
5 260 674
and what it writes by then: "nor its writer were in fact very fine a young man too like you whose very countenance may vouch for the honour of her situation in life is a question which mr darcy is not possible for anyone in particular places at his gallantry but there was such a comfort"
Run that five generations deep over eight independent chains and the numbers fall in step: the average vocabulary drops from 1,070 words to 260, and the count of distinct two-word states from 24,870 to 674. What comes out stays grammatical and loses its footing slowly — one of the fifth-generation chains writes: “nor its writer were in fact very fine a young man too like you whose very countenance may vouch for the honour of her situation in life”.
Figure 4: Five generations of a model trained on its own output: the vocabulary and the number of distinct two-word states both fall away.
Five generations is a small demonstration of a real effect. Ilia Shumailov and colleagues showed in Nature in 2024 that training models recursively on their own output degrades the distribution — the rare cases thin out first, then the language narrows toward what is common — and that keeping the original human data in the mix substantially reduces the damage.
The lesson is not that self-training is fatal. It is that a chain can only visit states it has been given, so a model trained on a narrowed world stays in that world.
Source: Claude E. Shannon, “A Mathematical Theory of Communication,” Bell System Technical Journal 27, 1948, 379–423 and 623–656; Ilia Shumailov et al., “AI models collapse when trained on recursively generated data,” Nature 631, 2024, 755–759.
Seven shuffles
The last stop is the one Ulam started from, because a shuffled deck is a chain too. Every arrangement of the cards is a state, and every shuffle is a step.
The shuffle worth modelling is the dovetail: cut the deck at a random point, then interleave the two halves so that the order inside each half is kept. Dave Bayer and Persi Diaconis worked out its mixing time in 1992, and got the number that has been repeated ever since — seven.
Seven is the point where the deck is close to random, not where it becomes random. Their table gives the total variation distance from a well-shuffled deck as 0.924 after five shuffles, 0.614 after six, 0.334 after seven and 0.167 after eight. The distance is the fraction of the deck’s arrangements that would have to be re-weighted to make the shuffle exactly uniform, so 0.334 is still a noticeable gap.
Figure 5: Distance between the shuffled and random distributions of rising sequences, by shuffle count.
You can see the same thing without their mathematics by watching one statistic: how many rising sequences the deck has. A deck in perfect order has one. A random deck has about 26.5, since the average is half of 52 plus a half. Each riffle can at most double the number.
"""How many riffle shuffles randomise a deck? Measured, not quoted."""import numpy as npN, TRIALS =52, 30000def riffle(deck, rng):"""Cut at a Binomial(52, 1/2) point, then interleave, keeping each packet's order.""" cut = rng.binomial(N, 0.5) left, right = deck[:cut], deck[cut:] out, i, j = np.empty(N, dtype=np.int8), 0, 0for k inrange(N):if rng.random() < (left.size - i) /max(1, (left.size - i) + (right.size - j)): out[k] = left[i]; i +=1else: out[k] = right[j]; j +=1return outdef rising(deck):returnint((np.diff(deck) <0).sum()) +1def spread(values): counts = np.bincount(values, minlength=N +1)[1:N +1]return counts / counts.sum()rng = np.random.default_rng(1992)base = spread(np.array([rising(rng.permutation(N)) for _ inrange(TRIALS)]))print(f"{'shuffles':>9}{'top card still on top':>23}{'rising sequences':>18}{'distance from random':>22}")for n inrange(1, 11): tops, rs =0, np.empty(TRIALS, dtype=np.int16)for trial inrange(TRIALS): deck = np.arange(N, dtype=np.int8)for _ inrange(n): deck = riffle(deck, rng) tops +=int(deck[0] ==0) rs[trial] = rising(deck) dist =0.5*float(np.abs(spread(rs) - base).sum())print(f"{n:>9}{tops / TRIALS:>23.2%}{rs.mean():>18.2f}{dist:>22.3f}")print(f"\na random deck keeps the same top card {1/52:.2%} of the time")
shuffles top card still on top rising sequences distance from random
1 50.58% 13.76 0.999
2 24.59% 20.12 0.882
3 12.54% 23.31 0.548
4 6.39% 24.91 0.297
5 3.85% 25.70 0.147
6 2.88% 26.09 0.074
7 2.34% 26.29 0.036
8 2.05% 26.41 0.016
9 1.91% 26.44 0.015
10 1.91% 26.47 0.010
a random deck keeps the same top card 1.92% of the time
The distance falls by roughly half with each shuffle — 0.999, 0.882, 0.548, 0.297, 0.147, 0.074, 0.036 — and by seven shuffles there is almost nothing left for this statistic to measure. The top card is back to being the top card about 2% of the time, which is the random rate.
Two honest caveats about that table. The distance reported is between the distributions of one statistic, not between the full distributions of deck arrangements, so it is a lower bound on the true distance — which is why it reaches 0.036 at seven shuffles while Bayer and Diaconis give 0.334. And my riffle is the tidy textbook one. A real riffle, done by hands, mixes faster or slower depending on how well the packets interleave, and the common estimate that an overhand shuffle needs thousands of passes comes from a different and much slower model.
Where the chain stops
Everything above works because the state is small and the rules do not change. Both of those are assumptions, and both fail somewhere useful.
The first failure is memory. A chain that conditions on the last two words cannot know what the sentence is about, and one that conditions on the current temperature cannot know that the planet has been warming for a century. Where the thing being modelled feeds back into itself — water vapour rising with temperature, models trained on text that models wrote — the state is not enough, and the chain gives a confidently narrow answer.
The second failure is stationarity, the assumption that the table of moves stays the same. A chain estimated on last decade’s data will keep predicting last decade’s behaviour, because nothing in it can represent a change in the rules.
Those are not reasons to abandon the idea. They are the reasons it is still being extended: hidden states for the things you cannot observe, time-varying tables for the things that drift, and attention for the history a single state cannot hold. The pattern is always the same. Find the state, write down the moves, and multiply.
Which is a good place to leave Markov’s feud. He won the argument he was having — independence is not required for regularity, and his two numbers from a poem prove it — and then the machinery he built to win it went further than the argument ever did.
Sources and further reading
The poem and the feud:
A. A. Markov, An Example of Statistical Investigation of the Text Eugene Onegin Concerning the Connection of Samples in Chains (1913), English translation in Science in Context 19(4), 2006, 591–600 — the translation; original in the Bulletin of the Imperial Academy of Sciences of St Petersburg 7(3), 153–162 — record
A. A. Markov, Extension of the Law of Large Numbers to Quantities Depending on Each Other (1906), Izvestiya of the Kazan Physico-Mathematical Society
Eugene Seneta, Markov and the Birth of Chain Dependence Theory, International Statistical Review 64, 1996, 255–263 — record
Eugene Seneta, Statistical Regularity and Free Will: L. A. J. Quetelet and P. A. Nekrasov, International Statistical Review 71, 2003, 319–334
Jacob Bernoulli, Ars Conjectandi, 1713
Adolphe Quetelet, A Treatise on Man and the Development of his Faculties (Edinburgh, 1842), the passage on the budget of prisons, dungeons and scaffolds, pp. 82–83 — scan; French original, Sur l’homme et le développement de ses facultés (Paris, 1835) — scan
