A restaurant and a shuffled deck

How many new kinds the next sample will show estimated, with no model at all, how many species remain to be found. Its method runs out soon, since the first sample says nothing about what it has not seen. Going further needs an assumption about how kinds are distributed, and the most natural such assumption has an extraordinary range. It gives the alleles of a gene under neutral evolution, it gives the species in a sample, and in one special case it gives the cycles of a shuffled deck. It gets the shares of the prime factors of a random number too, although nothing there is random in the same way.

Warren Ewens found the formula in 1972 while studying genetics. A little later a memorable picture was attached to it, by Lester Dubins and Jim Pitman, and it is the easiest way into the whole structure. Customers come into a restaurant with infinitely many tables. The first sits at a table. Each customer after that either joins an occupied table, with probability proportional to the number of people there, or starts a new one with probability proportional to a fixed number $\theta$. That is all there is to it.

Why the number of tables is so simple #

Customer $i + 1$ starts a new table with chance exactly $\theta/(\theta + i)$, whatever the arrangement of the earlier customers. That is the key. The denominator is the total weight of all the choices, $i$ people and $\theta$ for a new table, and it does not depend on how those $i$ people have divided themselves. So “a new table at step $i$” is a sequence of independent events. The number of tables is a sum of independent coin tosses with shrinking chances, and:

$$\mathbb{E}[K_n] = \sum_{i=0}^{n-1} \frac{\theta}{\theta + i} \approx \theta \ln\left(1 + \frac n\theta\right).$$

At 10,000 customers, the exact sum gives 5.59, 17.58 and 69.60, and simulations agree within their sampling error. It grows like $\theta \ln n$, so that a thousand times more data brings only about $7\theta$ more kinds. There is always another new one to be found, ever more rarely, with a chance of $\theta/(\theta + n)$ at the next draw. That is a clean answer to the earlier essay’s question within this model. The singleton count from Good and Turing’s method estimates the same thing, as the expected number of tables with exactly one customer is about $\theta$, whatever $n$ is, so the estimate of “next is new” is about $\theta/n$. Good and Turing’s formula does not assume any of this, and yet the two agree where the model holds.

Ronald Fisher’s logarithmic series of 1943, for butterflies collected by Corbet in Malaya, had already found the same shape: numbers of species represented by $j$ individuals fall like $\theta x^j/j$, and the total grows like a logarithm. Fisher’s $\alpha$ is this $\theta$. So ecologists had been using the parameter thirty years before Ewens, and for quite a different reason.

Twenty-three at a table #

The scheme has an unexpected link to the oldest question about collisions. All $n$ customers sit alone with probability $\prod_{i<n} \theta/(\theta + i)$, which for large $\theta$ is close to $\exp(-n^2/2\theta)$. That is the birthday calculation, with $\theta$ playing the part of the number of days. Put $\theta = 365$. The chance that 23 people all sit apart is 0.5070, and the same 23 in the usual model of 365 days give 0.4927, so the even-odds point moves only from 23 people to 24. They differ in that here a repeat is chosen in proportion to popularity rather than uniformly, and that matters less than one might think, since a first repeat is a matter of pairs. So the restaurant with large $\theta$ can be seen as a birthday problem in which the number of possible birthdays is infinite and the distribution comes from the sampling itself.

It also gives a simple way of reading $\theta$ from a sample, by setting the observed number of kinds equal to the sum. A catch of 1,000 individuals with 40 species gives $\theta = 8.21$, and it predicts that doubling the catch brings about 6 more species, whereas a million individuals would bring only about 57 more. That prediction is much stronger than anything the assumption-free extrapolation can make, and it is only as good as the model.

The cycles of a shuffle #

Now the surprise. Set $\theta = 1$ and compare the tables with the cycles of a random permutation.

The bars match, and they lie along $1/j$: the number of fixed points of a shuffle averages one, the number of 2-cycles averages one half, and so on. In 4,000 runs the cycle counts are 1.020, 0.488 and 0.330, and the tables give 1.015, 0.502 and 0.335. It is not an approximation. It is an exact identity, as the following construction shows. Build a permutation in cycle notation, one number at a time. Number $i + 1$ can either close off into a new cycle on its own, in one way, or be inserted after any of the $i$ earlier numbers, which puts it into that number’s cycle. So of $i + 1$ equally likely places, one starts a new cycle, and a cycle is chosen with probability proportional to its length. That is the restaurant at $\theta = 1$, and it is the Feller coupling, described by William Feller in 1945.

So everything known about cycles carries over. The fixed points of a hundred prisoners’ shuffle, and the fact that the longest cycle exceeds half with probability $\ln 2$, are statements about a restaurant. The number of cycles is the number of tables: about $\ln n$. It also generalises, because for $\theta \ne 1$ it is still a distribution on permutations: give each permutation a weight of $\theta$ to the power of its number of cycles. For an exact check, all 5,040 permutations of seven elements were grouped by their cycle type, and each share is the sampling formula at $\theta = 1$ to twelve places.

One formula for every split #

The probability of any particular arrangement, with $a_j$ tables of size $j$, is Ewens’s sampling formula:

$$P(a_1, \ldots, a_n) = \frac{n!}{\theta(\theta + 1)\cdots(\theta + n - 1)} \prod_{j} \frac{\theta^{a_j}}{j^{a_j} a_j!}.$$

At $\theta = 1$, six draws all at one table have probability 0.1667, one in six, and all at separate tables it is $1/720$. Larger $\theta$ favours more kinds. The eleven values add up to one for all $\theta$, and the sum was checked at sizes 5, 8 and 12 for three values of $\theta$. Look at how $\theta$ appears in the formula. Only through $\theta^K$, where $K$ is the number of tables. So given $K$, the probability of a split does not depend on $\theta$ at all, which means that the number of kinds is a sufficient statistic for $\theta$. Everything the sample has to say about $\theta$ is contained in how many kinds there are, and the rest of the configuration tells nothing.

That is what made the formula so useful in genetics. Motoo Kimura’s neutral theory holds that most variation between alleles has no effect on fitness. Given the number of alleles observed in a sample, the distribution of their counts is then completely known, and it does not involve the unknown mutation rate. So one can compare the actual counts against that exact distribution, as a test of neutrality. Ewens’s formula gives this directly. Mutations in this model are new tables, and $\theta = 4 N_e \mu$ measures the rate at which they come, in units of the population size.

John Kingman showed in 1978 that the scheme is, in a sense, the only one of its type: a partition of a sample is exchangeable and consistent as the sample grows exactly when it comes from a random “paintbox”, in which each individual picks a colour according to a fixed random list of shares. So there has to be a limit distribution for the shares, and it has a tidy description. Break a stick, keeping a $\mathrm{Beta}(1, \theta)$ fraction, then break what remains in the same way, and sort the pieces. That is the Poisson–Dirichlet distribution.

The largest table #

Given that there are about $\theta \ln n$ tables, a guess would be that any one table holds a vanishing share of everyone. It does not. The largest table keeps a positive share however large $n$ becomes, and the same goes for the second and third. At $\theta = 1$ the simulation gives 0.6265 for the largest, near the limit of 0.62433, the Golomb–Dickman constant. Lawrence Shepp and Stuart Lloyd in 1966 found this as the mean length of the longest cycle of a random permutation.

Then comes the most striking connection of all. Take a random number near $N$, and consider its prime factors. Karl Dickman studied in 1930 the size of the largest prime in a typical number, and that size, expressed as a power of $N$, has a limit distribution. Its mean is 0.62433, the same constant. Patrick Billingsley showed in 1972 that the whole list is the same: the logarithms of the prime factors of a random number, divided by $\ln N$ and sorted, converge to Poisson–Dirichlet with $\theta = 1$. So a random number’s factorisation looks statistically like a shuffled deck’s cycles. There is a real reason, as Richard Arratia, Andrew Barbour and Simon Tavaré explored. Both structures decompose an object uniquely into irreducible pieces, a cycle or a prime, and in both cases the number of irreducible objects of size $j$ is about $1/j$ of the available weight. Polynomials over finite fields behave in the same way, and the way they factor has precisely this limit, and the same is true of the components of random mappings, at $\theta = 1/2$. So the numbers follow the distribution of a card shuffle, even though nothing is being shuffled, and this has become a useful heuristic for making guesses about how many primes a typical number has.

A second parameter #

One thing the Ewens scheme cannot produce is a power law. Real collections of kinds frequently have far more than $\theta \ln n$ kinds, and they have heavy-tailed sizes. Words in English are the standard example, with their vocabulary growing like a power of the length of text.

Subtract a small discount $\sigma$ from every occupied table’s weight and add it to a new table’s weight. Then the new-table chance is $(\theta + \sigma K)/(\theta + n)$, which is no longer independent of the past. It grows as tables accumulate, and the feedback changes the growth. Pitman and Marc Yor’s two-parameter scheme gives $K_n \sim n^\sigma$ for $0 < \sigma < 1$. The straight lines on the figure’s logarithmic axes have roughly this slope, with 56, 338 and 2,769 kinds at 30,000 customers. The singletons are now a fixed share of everyone, at 0.0048 to 0.0696 of customers, as opposed to a constant number. That is the signature of a power law, and it is a relative of the digit laws that come from multiplying, in that a simple self-reinforcing mechanism gives a heavy tail.

A further feature of this family deserves mention, since it decides whether a sample could tell the two kinds of process apart. Ewens’s scheme has singletons in a fixed number whatever the sample size, while a discount has them in a fixed share. So a simple test is to double the sample and look at how many singletons there are. Doubling the catch above leaves its expected 8 or so singletons unchanged, under the Ewens model, though the sample has grown. Under a discount they would double. It is a very direct way of distinguishing the logarithm from the power, and it needs no fitted parameter of any kind at all, and a small survey can do it.

These schemes have been widely used in statistical language modelling, and they fit word counts remarkably well in the middle of the range. Nobody thinks that writers are sitting at tables, of course. The fit shows that the statistics of new words are close to those of a rich-get-richer process with a small new-word cost. This has been said many times about many things, and it is worth remembering that a good fit of a curve is only weak evidence for a mechanism.

How far a simulation can go #

The exact parts here are the sum for tables, the formula, its check on small sizes, and the identity for permutations. The rest is sampling. Shares and counts come from modest numbers of seeded runs, with sampling error visible in the last digits. The largest-table figure at 3,000 customers has also not fully converged to its limit, so that it could differ by a part in a hundred from the true value. The power-law slopes cannot be read exactly, since convergence is very slow: 30,000 customers is not a large sample for an exponent, with a value of 0.75 giving about 2,769 kinds where a pure power would give somewhat more, because the slope settles slowly.

Nothing here is a claim about a real ecology or genome. A sample that fits the formula has passed one test, and that is consistent with neutrality without demonstrating it. Equally, the numbers-and-permutations analogy is a theorem about limits, and it says nothing about any particular number.

What was computed and what was quoted #

The sum for the mean, and the formula’s total, were computed directly. All 5,040 permutations of seven were enumerated for the check at $\theta = 1$. The other figures are seeded simulations: restaurants by direct sampling, cycles by a Fisher–Yates shuffle, with both sizes given. The discount scheme uses a sampler that picks an existing customer’s table and thins it by rejection, and that is exact in distribution.

Quoted: the theorems of Ewens, Kingman, Pitman and Yor, Shepp and Lloyd, Dickman, Billingsley, and Arratia, Barbour and Tavaré; Fisher’s series; the genetics interpretation; and the constant 0.62433. The identity between restaurant and permutation is proved above, and the one between restaurant and prime factors is quoted only.

Still open: where the analogy with primes stops #

How far does the analogy between prime factors and cycles extend? At the level of the largest few shares it is a theorem. Finer questions are open. Do the prime factors of integers near a given size have the same joint statistics at a very fine scale as a random permutation’s cycles? Whether the two match in the small print, in questions about whether two numbers share a divisor in a given range, or the fine-scale distribution of the number of divisors, is a matter of active study by Kannan Soundararajan, Adam Harper, Andrew Granville and others. Harper showed that a certain random multiplicative model has behaviour that is better than expected, and it is not yet clear where the permutation analogy gives a wrong answer.

For word statistics the open question is different. What mechanism produces the observed exponents? Herbert Simon’s rich-get-richer model of 1955 and many successors reproduce the form, and none explains why language should have this particular exponent, and it is not clear that a single explanation exists.

One rule, and shares that stay #

A customer joins a table with probability proportional to its size or starts a new one with weight $\theta$, and everything follows from one exact formula: about $\theta \ln n$ tables, a number of kinds that is sufficient for $\theta$, and shares that stay positive however many customers come. At $\theta = 1$ the restaurant is exactly the cycle structure of a shuffled deck.

And the same shares then appear in the prime factorisation of a random number. The largest table, the longest cycle and the largest prime take about 62% of the whole on average, on their respective scales. It is not common for one constant to tie together something as concrete as a shuffle, something as biological as a neutral gene and something as rigid as an integer’s factors. It seems to be just what happens when a large object breaks uniquely into irreducible pieces, and nobody has a theorem saying how general that is.