charlesreid1.com blog

Baby Mac Crypto

Posted in Mathematics

permalink

The fictional Acme Insurance Corp employs a team of recreational mathematicians, people who study the fun parts of math (game theory, combinatorics, and cryptography), and puts them to work as actuaries. They take the most entertaining problems in mathematics and turn them into dry, boring actuarial tables.

This week, two fast food chains want Acme to insure the prizes in their collect-the-set contests. Burger Slop and Chicken Star will each print ten million game pieces, and almost none of them will win anything big. When a customer completes a set, Acme pays the prize. Each chain is choosing among three versions of its game, so Acme needs to know how much it should expect to pay out for each version, and which version it would rather insure.

Logos of the two fictional chains: Burger Slop, a brown hexagon lying on a flat side, and Chicken Star, a navy star with rounded corners, each with its name in white

Throughout, we assume that no prize is ever stolen. The theft section near the end asks what happens when that assumption fails.

The problem is small enough to solve on paper, and solving it turns up three nuggets:

  1. Collecting a set is a coupon collector problem where one coupon is tens of thousands of times rarer than the others.
  2. The player holding a rare piece is not a typical player, because a player with more pieces is more likely to hold it. Treating them as typical can underprice a game by a factor of about three.
  3. All the information needed to check a code is contained in the code itself.

The Contests

Burger Slop's and Chicken Star's games work the same way. Each game piece shows one property on a game board, printed in its set's color, and a 14-character code. The properties are grouped into sets, and a player who holds a piece of every property in a set wins that set's prize by entering the pieces' codes online. Each code counts once, for the first player to enter it. Each set has one rare property and some common ones. "Common" here means plentiful: each common property is printed hundreds of thousands of times, while only as many pieces carry the rare property as the set has prizes.

All three versions of the game use a board of 24 properties, split into sets of \(k\) properties: one rare and \(k-1\) common.

game \(k\) sets
duopoly 2 12
triopoly 3 8
tetropoly 4 6

Three boards of 24 squares: duopoly cut into 12 sets of 2, triopoly into 8 sets of 3, tetropoly into 6 sets of 4, each set in its own color, with one rare property starred in each set

Each chain prints \(N = 10^7\) pieces. Each prize is one rare piece: 200 prizes of $1,000, 20 of $25,000 and one grand prize of $1,000,000, for $1.7 million in all. Common pieces make up 60% of the pieces, and the rest carry no property.

The board's 24 properties make \(24/k\) sets, each with \(k-1\) common properties: 12 common properties in duopoly (12 sets of 1), 16 in triopoly (8 sets of 2) and 18 in tetropoly (6 sets of 3). They split the 60% equally, so each common property has frequency

$$ p = \frac{0.6}{(24/k)(k-1)} $$

which is 1 in 20 for duopoly, 1 in 26.7 for triopoly and 1 in 30 for tetropoly.

The contests run six weeks, and half the pieces go to people who never enter a code. The rest go to players, and here the chains differ. At Burger Slop, 80% of players buy little and collect 3 pieces over the contest, and the other 20%, the heavy players, buy a lot and collect 30. At Chicken Star, every player collects one piece a week, so 6 over the six-week contest. That makes about 600,000 players at Burger Slop and 830,000 at Chicken Star. Each piece is an independent random draw from the printed pieces, players do not trade, and finding a rare piece does not change how many more pieces a player collects.

A Skewed Coupon Collector

In the classical coupon collector problem, each draw is uniform over \(k\) coupons, and collecting all \(k\) takes \(k H_k\) draws on average, where \(H_k = 1 + \tfrac{1}{2} + \cdots + \tfrac{1}{k}\) is the \(k\)-th harmonic number, which grows like \(\ln k\). These coupons are about as far from uniform as coupons get.

Model each piece a player receives as an independent draw from the \(N\) printed pieces. A set \(S\) has one rare piece and \(k-1\) common pieces. Let \(m_S\) be the number of prizes for \(S\). The rare piece is printed \(m_S\) times, so a draw is the rare piece with probability

$$ \epsilon_S = \frac{m_S}{N} $$

There are 221 prizes in all. Even a set that held every one of them would have \(\epsilon_S = 221/10^7\), so its rare piece would turn up once in \(10^7/221 = 45{,}249\) draws on average. A real set has fewer prizes, so its rare piece is scarcer and the wait is longer. This is the same in every game, because it depends only on how often the rare piece is printed.

Let \(T_S\) be the number of draws until a player holds all \(k\) pieces of \(S\), \(T_R\) the number until the rare piece, and \(T_C\) the number until all \(k-1\) common pieces. A player holds the set once they hold both, so \(T_S = \max(T_R, T_C)\) and \(T_R \le T_S \le T_R + T_C\). Averaging,

$$ \frac{1}{\epsilon_S} \le E[T_S] \le \frac{1}{\epsilon_S} + E[T_C] $$

The last term is the average wait for the common pieces, and this is where the games differ:

  • Duopoly: one common piece per set, turning up 1 draw in 20. It takes 20 draws on average.
  • Triopoly: two common pieces, each 1 in 26.7. Either one will do for the first, so the first turns up twice as fast, after \(26.7/2 = 13.3\) draws. The second takes another 26.7, for 40 in all.
  • Tetropoly: three common pieces, each 1 in 30. The first takes \(30/3 = 10\) draws, the second another \(30/2 = 15\), and the last another 30, for 55 in all.

Each total is \(E[T_C] = H_{k-1} \times (1/p)\), where \(1/p\) is the wait for one particular common piece: \(1 \times 20 = 20\), \(1.5 \times 26.7 = 40\) and \(1.83 \times 30 = 55\). So the common pieces add 20 to 55 draws to a wait of at least 45,249 for the rare piece. The rare piece's wait is more than 800 times longer, in every game. The 20 to 55 draws are an upper limit, because players collect the common pieces while they wait for the rare one and almost always have them all by the time it shows up.

So the wait for any set is at least 45,249 draws, and all of it is spent waiting for the rare piece. No player waits that long, though: players collect 3 to 30 pieces. What the rest of the post uses is the same fact for a fixed number of draws. Draws are independent, so finding the rare piece says nothing about a player's other pieces, and the player completes the set exactly when those include every common property.

One Number

A prize is paid only if its rare piece is redeemed, so everything the actuary needs comes down to one number: \(q\), the probability that a given rare piece is redeemed. That takes two things. The piece has to go to a player, which happens with probability 1/2. And that player's other pieces have to include all \(k-1\) common properties of its set.

Size Bias

In the figure, each person is 10% of a chain's players and each square is one game piece over the six-week contest.

Under the Burger Slop logo, players in two groups: eight people with 3 pieces each, and two heavy players with 30 pieces each. Under the Chicken Star logo, players in one group: ten people with 6 pieces each

The naive move is to take a typical Burger Slop player and ask how often their other pieces complete the set. But the rare piece does not go to a typical player. Every piece is equally likely to be the rare one, so a player with 30 pieces is ten times as likely to hold it as a player with 3. This is the inspection paradox.

Let \(X\) be the number of pieces a random player collects, so \(\Pr(X = n)\) is the fraction of players who collect \(n\) pieces. At Burger Slop, \(\Pr(X = 3) = 0.8\) and \(\Pr(X = 30) = 0.2\). At Chicken Star, \(\Pr(X = 6) = 1\). Weighting each player by their number of pieces, the holder of a rare piece received \(n\) pieces with probability

$$ \frac{n \Pr(X = n)}{E[X]} $$

The factor \(n\) is the weighting, and dividing by \(E[X]\), the average number of pieces per player, makes the probabilities add up to 1. How much this changes the answer depends on how unevenly the pieces are spread among a chain's players, and the two chains are opposites.

At Burger Slop the spread is wide. The average player collects

$$ E[X] = 0.8 \cdot 3 + 0.2 \cdot 30 = 8.4 $$

pieces, and the holder of a rare piece is a heavy player with probability \(0.2 \cdot 30 / 8.4 = 0.714\). Heavy players are 20% of players but hold 71% of the rare pieces. Averaging \(n\) with these weights, the holder of a rare piece has \(E[X^2]/E[X] = 22.3\) pieces, against 8.4 for a typical player.

At Chicken Star there is no spread. Every player collects 6 pieces, so \(E[X] = 6\), and the holder received 6 pieces with probability \(6 \cdot 1 / 6 = 1\). The holder of a rare piece is a typical player, and the naive answer is exact.

Collecting the Commons

The holder of a rare piece wins only if their other pieces include every common property of its set. Write \(P_k(r)\) for the chance that \(r\) draws include all \(k-1\) common properties of a set, each of frequency \(p\).

In duopoly there is one common property. A draw misses it with probability \(1-p\), so \(r\) draws all miss it with probability \((1-p)^r\), and

$$ P_2(r) = 1 - (1-p)^r $$

With two or more common properties, subtracting the chance of missing each one subtracts too much: draws that miss two properties at once get subtracted twice. Inclusion-exclusion, the same tool as in Enumeration Solitaire, adds those back, then corrects for draws that miss three, and so on. A run of \(r\) draws misses a given group of \(a\) properties with probability \((1-ap)^r\), so for triopoly and tetropoly

$$ P_3(r) = 1 - 2(1-p)^r + (1-2p)^r $$
$$ P_4(r) = 1 - 3(1-p)^r + 3(1-2p)^r - (1-3p)^r $$

Putting the pieces together gives the one number. A rare piece reaches a player with probability 1/2, its holder collected \(n\) pieces with the size-bias weight \(n \Pr(X = n) / E[X]\), and the holder's other \(n - 1\) pieces include every common property with probability \(P_k(n-1)\). Summing over \(n\) gives

$$ q = \frac{1}{2} \, \frac{E\left[ X \, P_k(X-1) \right]}{E[X]} $$

Now put in Burger Slop's numbers. A heavy player collects 30 pieces, so a heavy player holding a rare piece has 29 others. They include the commons the player needs with probability \(P_k(29)\): 0.77 in duopoly, 0.44 in triopoly and 0.24 in tetropoly. A light player collects 3 pieces, so a light player holding a rare piece has only 2 others. They include the commons with probability \(P_k(2)\): 0.10 in duopoly, under 0.003 in triopoly, and never in tetropoly, which needs three.

With the weights 0.286 and 0.714 from the size bias section, the formula for Burger Slop is

$$ q = \tfrac{1}{2} \left( 0.286 \, P_k(2) + 0.714 \, P_k(29) \right) $$
game \(q\) \(q\) without size bias \(q\) simulated
duopoly 0.290 0.116 0.293 ± 0.004
triopoly 0.159 0.046 0.156 ± 0.004
tetropoly 0.084 0.024 0.082 ± 0.003

The middle column uses a typical player's weights, 0.8 and 0.2, in place of 0.286 and 0.714. The simulation deals pieces to players and counts the rare pieces that end up in completed sets. It agrees with the formula, and not with the middle column. Ignoring size bias underprices every Burger Slop game by a factor of 2.5 to 3.5.

Bigger sets lower \(q\) twice over: the holder must find more common properties, and each is scarcer, because more of them share the same 60% of the pieces.

Chicken Star

At Chicken Star the holder of a rare piece has 5 other pieces, so

$$ q = \tfrac{1}{2} P_k(5) $$
game Burger Slop Chicken Star
duopoly 0.290 0.113
triopoly 0.159 0.013
tetropoly 0.084 0.001

Chicken Star's \(q\) is lower in every game, and far lower in the bigger ones: with 5 other pieces, a holder rarely finds all three common properties of a tetropoly set.

Contest Length

The formula for \(q\) never mentions how long the contest runs. A set counts if its holder has every piece by the end, so what matters is how many pieces a player collects, not when they collect them.

At Burger Slop, say light players collect 0.5 pieces a week and heavy players 5. Over \(w\) weeks they collect \(0.5w\) and \(5w\) pieces, whole numbers for contests of 4, 6 and 8 weeks, and

$$ q = \tfrac{1}{2} \left( 0.286 \, P_k(0.5w - 1) + 0.714 \, P_k(5w - 1) \right) $$

The weights are still 0.286 and 0.714, because both rates scale with \(w\). At six weeks this is the \(q\) from before, with 3 and 30 pieces.

Whether the length cancels out depends on how the chain prints. If it prints a fixed run of \(N\) pieces and runs the contest until they are gone, it does: \(N\) and the number of customers set the pieces per player, and \(w\) drops out. If it prints for the length of the contest, it does not: at the same buying rates, players collect more pieces in a longer contest, with the same 221 prizes.

Here is Burger Slop's \(q\) at four, six and eight weeks:

game 4 weeks 6 weeks 8 weeks ratio, 8 to 4
duopoly 0.230 0.290 0.329 1.4
triopoly 0.093 0.159 0.214 2.3
tetropoly 0.035 0.084 0.138 3.9

Stretching the contest from four weeks to eight multiplies duopoly's \(q\) by 1.4 and tetropoly's by 3.9.

At Chicken Star every player collects one piece a week, so over a contest of \(w\) weeks each player collects exactly \(w\) pieces, and

$$ q = \tfrac{1}{2} P_k(w - 1) $$
game 4 weeks 6 weeks 8 weeks ratio, 8 to 4
duopoly 0.071 0.113 0.151 2.1
triopoly 0.004 0.013 0.024 6.0
tetropoly 0.0001 0.0010 0.0032 29

Chicken Star's \(q\) is far more sensitive to the length than Burger Slop's. Its players hold too few pieces to finish a big set, so each extra week raises their chances by a large factor: from four weeks to eight, tetropoly's \(q\) grows 29-fold.

Expected Payout

If Acme Insurance Corp takes on a contest, it pays every prize that is won. Each prize is won with probability \(q\), so the prize money Acme expects to pay out is \(q\) times the $1.7 million at stake. The premium Acme would charge the chain starts from this expected payout and adds a margin for overhead and profit; this post stops at the expected payout.

Because \(q\) depends on how long the contest runs, so does the expected payout. Here it is for contests of four, six and eight weeks.

At Burger Slop:

game \(q\), 4 weeks payout, 4 weeks \(q\), 6 weeks payout, 6 weeks \(q\), 8 weeks payout, 8 weeks
duopoly 0.230 $390k 0.290 $494k 0.329 $560k
triopoly 0.093 $158k 0.159 $270k 0.214 $365k
tetropoly 0.035 $60k 0.084 $143k 0.138 $235k

At Chicken Star:

game \(q\), 4 weeks payout, 4 weeks \(q\), 6 weeks payout, 6 weeks \(q\), 8 weeks payout, 8 weeks
duopoly 0.071 $121k 0.113 $192k 0.151 $256k
triopoly 0.004 $6.9k 0.013 $21k 0.024 $42k
tetropoly 0.0001 $190 0.0010 $1.7k 0.0032 $5.4k

The expected payout holds only for a contest of a stated length, so Acme's policy should state the length.

The expected payout is not what Acme will actually pay. Two rare pieces almost never share a holder, so the prizes are paid nearly independently, and the grand prize is most of the money. The payout is mostly one coin flip.

At Burger Slop, the 220 smaller prizes are expected to cost $203k in duopoly. So duopoly pays about $203k with probability 0.71, and about $1.2M with probability 0.29, when the grand prize is won too. Tetropoly's smaller prizes are expected to cost $59k, so it pays about $59k with probability 0.916 and about $1.06M with probability 0.084. Chicken Star's tetropoly pays the grand prize with probability 0.001, and its smaller prizes are expected to cost $700 in all.

Sensitivity

Every \(q\) above rests on a guess about how many pieces players collect. What matters to Acme is how far a wrong guess would leave its expected payout too low. So that the guesses can be compared, each one is moved by the same amount, a third.

At Burger Slop the main guess is the figure of 30 pieces per heavy player. Here is \(q\) if the true figure is a third lower or higher:

game \(q\) at 20 \(q\) at 30 \(q\) at 40 ratio, 40 to 20
duopoly 0.213 0.290 0.344 1.6
triopoly 0.082 0.159 0.230 2.8
tetropoly 0.031 0.084 0.149 4.8

Bigger sets are more sensitive. Let \(\lambda = pr\), the expected number of copies of one common property in \(r\) draws. Replacing each \((1-ap)^r\) in \(P_k(r)\) by \(e^{-a\lambda}\) gives \(P_k(r) \approx (1 - e^{-\lambda})^{k-1}\), so a 1% error in \(r\) moves \(P_k\) by about

$$ (k-1) \, \frac{\lambda e^{-\lambda}}{1 - e^{-\lambda}} $$

percent. The fraction falls from 1 toward 0 as \(\lambda\) grows, since extra draws help less once a player probably holds each common property, and the factor \(k-1\) counts the common properties that must all turn up. For a heavy player's 29 other pieces this is 0.44% in duopoly, 1.1% in triopoly and 1.8% in tetropoly.

In dollars, here is how much each expected payout would be too low if a guess is a third too low. A third more pieces for every player is the eight-week contest from the contest length section.

a third too low chain duopoly triopoly tetropoly
heavy players' pieces (40, not 30) Burger Slop $91k $121k $110k
every player's pieces (eight weeks, not six) Burger Slop $66k $94k $92k
every player's pieces (8, not 6) Chicken Star $64k $20k $3.7k

Actuary Recommendation

Acme prefers the game whose expected payout is off by the fewest dollars when a guess about the players is a third too low. The coin flip does not enter: across the many contests Acme insures, the flips average out, but a wrong guess about players is wrong in every contest priced with it.

With no theft, insure duopoly at Burger Slop and tetropoly at Chicken Star.

At Burger Slop, tetropoly has the lowest expected payout, $143k against duopoly's $494k, but the heavy players make the bigger games' prices hard to pin down. If players collect a third more pieces than guessed, Acme's expected payout for duopoly is too low by $66k to $91k, and for the bigger games by $92k to $121k. Duopoly costs the most, and its cost is known best.

At Chicken Star, with 6 pieces each, almost nobody completes a tetropoly set. Its expected payout is $1,700, and 8 pieces per player would add $3,700 to it. Tetropoly is both the cheapest game and the one whose cost is known best.

Codes That Check Themselves

Strip away the paper and the contest is codes and an oracle: each piece's 14-character code is a token, and the server is an oracle that says whether a code is valid and what it wins. Forgery means making up a valid code. Theft means taking one that already exists. Acme insures a contest only if the chain uses Acme's codes, which make forgery hopeless and can be checked without a stored list of valid codes.

A code is 14 characters over 32 symbols (A-Z and 0-9, minus the look-alikes 0, O, 1 and I), so it is a 70-bit number, an element of \(D = \{0, 1, \dots, 2^{70}-1\}\).

The chain's market is split into 100 zones of equal population. Zone \(z\) gets \(N_z = 10^5\) pieces, and its piece \(i\), for \(i = 0, 1, \dots, N_z - 1\), gets the index \(j = z N_z + i\), so the indices run from 0 to \(N - 1\). A small table per zone says which property each \(i\) carries, with the zone's rare pieces at the lowest \(i\).

Acme holds one secret key, and \(\pi\) is a published format-preserving cipher on \(D\) under that key: a permutation of \(D\) that looks random to anyone without the key. Format-preserving means it turns a 70-bit number into another 70-bit number, so the result still fits in 14 characters. Piece \(j\) gets the code \(\pi(j)\).

To check a code \(y\), the server decrypts it to \(j = \pi^{-1}(y)\). If \(j < N\), the code is piece \(i = j - z N_z\) of zone \(z = \lfloor j / N_z \rfloor\); the server reads its property from that zone's table and marks \(j\) as redeemed. Otherwise the code is rejected.

Forgery. Without the key, a guessed code decrypts to an essentially uniform element of \(D\). It passes as a printed piece with probability

$$ \frac{N}{2^{70}} \approx 8.5 \times 10^{-15} $$

and as a rare piece with probability \(221 / 2^{70} \approx 1.9 \times 10^{-19}\). The index uses \(\log_2 N \approx 23.3\) of the code's 70 bits, and the other 46.7 bits act as an authentication tag. Forgery is not worth pricing.

The codes do nothing against theft, which the theft section takes up.

Theft

Every expected payout above assumes nobody steals a prize. A thief always redeems what they steal, so for a stolen piece \(q\) becomes 1, as it does for a holder who keeps buying until they have the commons. A stolen grand prize costs the full $1 million instead of \(q\) times $1 million, so it adds \((1 - q)\) times $1 million to the expected payout:

game Burger Slop Chicken Star
duopoly $710k $887k
triopoly $841k $987k
tetropoly $916k $999k

Duopoly's expected payout misses by the least. At Burger Slop this agrees with the answer without theft. At Chicken Star, let \(\theta\) be the probability that the grand prize is stolen. With 8 pieces per player, not 6, duopoly's expected payout is too low by $64k + \(\theta\) × $887k and tetropoly's by $3.7k + \(\theta\) × $999k, so duopoly is the better game once \(\theta\) is above 0.54.

The codes cannot prevent a theft; they only stop forgery. Preventing theft takes physical and digital security: chain of custody for the printed pieces, and protection for the oracle and its key, since a leaked key exposes every rare code. After a theft, the zone a code decrypts to can flag a suspicious claim, such as a rare code from one zone redeemed with common codes from another.

Summary

When one coupon is tens of thousands of times rarer than the rest, the expected wait for the whole set is the wait for that coupon, and holding the others does not make it any likelier to turn up.

The holder of a rare piece is not a typical player, because choosing a person through something they hold weights each person by how much they hold. This is the inspection paradox. At Burger Slop it is worth a factor of about three in price. At Chicken Star, where every player holds the same number of pieces, it is worth nothing.

Encrypting the piece index with a format-preserving cipher makes each code its own proof of validity. The server keeps one key and small tables instead of a list of \(10^7\) codes, and the spare bits of the code bound the forgery rate.

Code

https://git.charlesreid1.com/charlesreid1/baby-mac-crypto/src/branch/main/calculations.py

References

Fermi Problems Are Quadrature Problems

Posted in Mathematics

permalink

Here's a Fermi problem you'll run into sooner or later, usually in a job interview, sometimes over drinks with someone who wants to see how you think:

How many chocolate bars are eaten each year in the United States?

The provenance goes back to Enrico Fermi's Chicago physics classes, where he'd ask students to estimate the number of piano tuners in the city, or the width of a nail's head in miles, or anything else where the point wasn't the number but the decomposition. That habit of mind has since been laundered through McKinsey and Google into the standard interview format, but the underlying move is Fermi's: break a hard question into smaller questions you can each guess within a factor of two, and let the errors mostly cancel.

Does Size Matter?

Before we start estimating, there's an elephant in the room that deserves to be named and then, deliberately, ignored. What counts as "a chocolate bar"? A fun-size Snickers you get in a Halloween bucket? A king-size Twix? A square of a Hershey's bar? A Costco-sized Toblerone? The answer changes the final count by an arbitrary factor of 1/2 or 2 or worse - it's the single variable this whole problem is most sensitive to, and it's not really a quantitative question at all. It's a linguistic knob. Turning it doesn't teach us anything about populations or quadrature; it just relabels what we're counting.

So we're going to do what physicists do when a problem has a knob like this: we're going to assume a spherical cow. In the spirit of Fermi, "eating a chocolate bar" is an idealized, dimensionless event - a discrete unit of chocolate consumption, roughly the size of whatever the reader pictures when they hear the phrase. We're not going to quibble about grams or servings or whether a Reese's cup counts. If your definition differs from mine by a factor of two, your final answer differs from mine by a factor of two, and that's fine. The interesting structure of this problem lives in the population and its interactions, not in the definition of the counting unit, and the machinery we're about to build works the same way whichever definition you pick.

With that noted and set aside, on to the actual estimation.

The Naive Approach and What It's Secretly Assuming

The textbook approach to the chocolate question looks like this. Guess that the average American eats maybe one chocolate bar every few days - call it 0.3 per day. Multiply by the population, roughly 340 million. Multiply by 365. You get something on the order of 40 billion chocolate bars per year.

That's your answer, and if you've picked reasonable numbers it's probably within a factor of two or three of the truth.

But look at what that multiplication is quietly assuming: one average American, one average day, one average rate. Every source of variation in the real population has been smoothed into a single number. If you push on the estimate, the assumption starts to feel thin. A five-year-old on Halloween is not eating the same amount of chocolate as an fifty-year-old on Halloween. Neither of them is eating the same amount on Halloween as they are in March. And it gets worse when you notice that the seasonal spike itself depends on who you are: kids own Halloween, adults own Valentine's Day, families with young children own Easter. (We can pretty safely assume that seniors will be flat across all of it.)

You can try to refine the estimate by splitting into age buckets, or by splitting into seasonal buckets. Either helps a little. But do them separately and you still miss the thing that actually matters, which is that age and season interact. A "seasonal multiplier" averaged over the whole population says everyone eats three times more chocolate on Halloween. That's false in a specific and important way: kids eat ten times more, adults eat only slightly more, and the average is a fiction that lives in neither group.

The interaction is real, and no amount of separately-refined marginal averages will recover it. We need a formalism that lets us write the joint structure down directly.

Population as a Joint Density (or, A Little Quadrature Never Hurt Anyone)

Here's the reframe. Think of the population as a density over some space of coordinates. Each person has internal coordinates \(\xi\) - age, income, dietary preferences, whatever matters for the problem - and lives in the external coordinate of time \(t\). Let \(n(\xi, t)\) be the number density of people: \(n(\xi, t)\, d\xi\) is how many people have coordinates in a little box around \(\xi\) at time \(t\). Let \(r(\xi, t)\) be the per-capita consumption rate: chocolate bars per person per unit time, for a person with coordinates \(\xi\) at time \(t\).

Then the total number of chocolate bars eaten in a year is just the integral of the product:

$$ \text{CBE} = \int_0^{1\,\text{yr}} \int_\xi r(\xi, t) \, n(\xi, t) \, d\xi \, dt $$

Two objects, cleanly separated. \(n\) says who exists. \(r\) says what they do. The product \(r \cdot n\) is chocolate bars per unit coordinate per unit time, and integrating it over the whole domain gives the total. (This is the population balance equation, a workhorse in chemical engineering, particle dynamics, and demography. We're borrowing the machinery, not reinventing it.)

Now the key observation: any numerical evaluation of that integral is a quadrature. Pick a grid of abscissae \((\xi_i, t_j)\), assign each grid cell a weight \(w_{ij}\) (in the simplest case, just the cell's width), and sum:

$$ \text{CBE} \approx \sum_{i,j} w_{ij} \, r(\xi_i, t_j) \, n(\xi_i, t_j) $$

That's the whole game. The naive Fermi estimate we started with - rate × population × 365 - is exactly this sum with a single term: one \(\xi\)-bin covering everyone, one \(t\)-bin covering the whole year, one average rate. It's a one-point quadrature over the joint density. Refining a Fermi estimate is refining a quadrature grid. Everything from the back-of-envelope guess up to a full numerical integration lives on the same continuum.

And crucially, this formulation does not require \(r\) to factor as an age-dependent piece times a seasonal piece. You fill in each grid cell independently, so any interaction between coordinates is preserved by construction. That's exactly the structure the marginal-average approach destroys.

There's one more idea worth borrowing from numerical analysis here. The whole art of Gaussian quadrature is that abscissae should not be evenly spaced - they should cluster where the integrand varies fastest. The same instinct applies to Fermi estimation: a good quadrature grid matches the shape of what you're integrating, not the shape of the axis you're integrating over. For the chocolate problem, the "time" axis has twelve months, but the integrand has three spikes: Halloween, Valentine's, Easter. Using twelve evenly-spaced monthly bins would spend most of your resolution on quiet months where the integrand is flat and boring. Concentrating your grid on the three candy-heavy windows plus one bin for "everything else" captures almost all the structure with a quarter of the work. That reduction - twelve bins to four - is not a shortcut; it's the correct grid for this integrand.

There's a subtlety worth naming here, though. When we collapse twelve evenly-spaced months into four unevenly-sized "seasons," we're doing something a little sneaky: we're letting the content of the year (which holiday, which age group's spike) reshape the time axis itself. What was a pure external dimension - clock time, marching uniformly forward - is being warped into an internal dimension that carries information about the integrand. That works cleanly here because there's no other process on the time axis we care about. As problems gain dimensions, though, that alignment can break.

(As a more involved example, if we were trying to approximate number of chocolate bars eaten by a population, but under conditions of continuous changes in population or demographics, then integrating the total population over time would need a properly-resolved time axis - lumping 46 weeks into one bin would be throwing away real information. The seasonal grid works for our problem, because the integrand's structure and external axis's structure happen to coincide. If they don't, keep them separate, and "pay" for the extra abscissas (with a bit more bookkeeping).

Working the Chocolate Problem With Interactions

With the framework in hand, the estimate becomes a table. Rows are age buckets, columns are seasonal windows, each cell holds the number of bars eaten by that group in that window.

For the age axis, four buckets are enough: kids (0-12), teens (13-19), adults (20-64), and seniors (65+). Population is roughly uniform over these bands, so with 340 million Americans we can allocate roughly 55M kids, 30M teens, 190M adults, and 65M seniors. (Uniform-in-age is a lie, but a small enough one that the drama of this problem lives in \(r\), not in \(n\). That itself is a useful diagnostic - it tells us where refining would and wouldn't help.)

For the time axis, four bins: a Halloween window (~2 weeks around Oct 31), a Valentine's window (~2 weeks around Feb 14), an Easter window (~2 weeks around the spring holiday), and the remaining ~46 weeks of the year lumped into "rest of year." Numbers below are bars per person per day, and we'll multiply through by bucket population and bin length at the end.

Halloween (14d) Valentine's (14d) Easter (14d) Rest (312d)
Kids (55M) 3.0 0.4 1.5 0.3
Teens (30M) 2.0 0.6 0.4 0.4
Adults (190M) 0.5 1.5 0.3 0.25
Seniors (65M) 0.2 0.3 0.2 0.15

Look at that matrix for a second. The kids row peaks on Halloween. The adults row peaks on Valentine's. Easter has a bump for kids and almost nothing for anyone else. Seniors are flat and low. No product of a row-vector and a column-vector produces this pattern - the matrix is not rank-1, and any factorization \(r(\xi, t) = r_{\text{age}}(\xi) \cdot s(t)\) would flatten these ridges into a smooth surface that gets every cell wrong. This is precisely the structure the joint formulation preserves and the marginal-refinement approach loses.

Summing (population × rate × days) cell by cell:

$$ \text{CBE} \approx \sum_{i,j} n_i \cdot r_{ij} \cdot w_j \approx 44 \text{ billion bars/year} $$

Which is, satisfyingly, in the same order of magnitude as the naive estimate from the first section (40 billion for the naive approach, 44 billion for the quadrature estimate). That's the usual outcome, and it's not a knock on the framework - it's a reminder that averaging over correlated variables often lands close to the truth by luck. What the framework buys you isn't necessarily a better number; it's a legible number. Every approximation is named and located in a specific cell. If you wanted to sharpen the estimate, you'd know exactly where to add resolution - split Halloween into "trick-or-treat night" and "the week after," split kids into "young enough to trick-or-treat" and "too cool for it" - and you'd be adding grid points to the ridges, which is exactly what a higher-order quadrature scheme does automatically.

That's the transferable idea, and it generalizes beyond chocolate. Any Fermi problem about a population - how many haircuts per year, how many gallons of coffee consumed per week, how many miles driven per day - is an integral of a per-capita rate against a number density over some joint coordinate space. The reason the standard "multiply the averages" trick works at all is that most people's mental \(r\) is smooth enough that a one-point quadrature suffices. When it isn't - when the coordinates interact, when the ridges matter - the framework tells you exactly where to add resolution and why. You stop guessing a single average and start choosing a grid.

References

Tags:    mathematics    fermi problems    quadrature    numerical methods    chocolate   

Enumeration Solitaire: Derangements with Four Copies of Each Rank

Posted in Mathematics

permalink

Enumeration Solitaire is a one-player card game with a two-sentence rule set and a genuinely nontrivial answer. You shuffle a deck, deal cards face up, and count "ace, two, three, ..., king" over and over as you go. If the rank you speak ever matches the rank of the card you just dealt, you lose. Get through all 52 cards without a match and you win.

It is a four-copy generalization of the classical probleme des rencontres (the derangement problem), and the solution is one inclusion-exclusion argument compressed into one generating function.

Our working notes: Enumeration Solitaire.

The Problem

This one comes from the same FiveThirtyEight Riddler column as Baby Poker (archived on the Wayback Machine): Can You Deal With These Card Game Puzzles?.

The rules:

  • Take a standard shuffled 52-card deck.
  • Deal cards face up, one at a time.
  • As each card is dealt, speak the next rank in the cycle: ace, two, three, ..., king, then restart at ace.
  • If the spoken rank ever matches the rank of the dealt card, you lose.
  • If all 52 cards are dealt without a match, you win.

Two questions:

  • What is the probability of winning a single game?
  • On average, how many games must you play before winning once?

Setting Up the Problem

Label the deal positions \(i \in \{0, 1, \dots, 51\}\). At position \(i\) the spoken rank is \(s(i) = i \bmod 13\). Let \(A_i\) be the event "the card dealt at position \(i\) has rank \(s(i)\)." A win is the event that none of the \(A_i\) occur.

The naive move is to note that \(P(A_i) = 4/52 = 1/13\) for every \(i\) and treat the 52 events as independent. That gives \((12/13)^{52} \approx 0.0155\), and by linearity of expectation the average number of matches is \(52 \cdot 1/13 = 4\), so a Poisson approximation gives \(e^{-4} \approx 0.0183\). These are two different wrong answers because the events are not independent: the four cards of each rank compete for the same four labeled slots.

The right tool for "probability that none of a bunch of correlated events occurs" is inclusion-exclusion.

Inclusion-Exclusion

$$ P(\text{win}) = P\!\left(\bigcap_i A_i^c\right) = \sum_{k=0}^{52} (-1)^k S_k, \qquad S_k = \sum_{|I|=k} P\!\left(\bigcap_{i \in I} A_i\right). $$

Model the deal as a uniformly random permutation of the 52 distinct cards. For a fixed size-\(k\) subset \(I\) of positions, group its positions by the rank they demand: let \(m_r\) be the number of positions in \(I\) with \(s(i) = r\). Then \(m_0 + \dots + m_{12} = k\) and each \(m_r \le 4\). Call \((m_0, \dots, m_{12})\) the multiplicity profile of \(I\).

To build an arrangement in which every position of \(I\) is a match:

  1. For each rank \(r\), fill the \(m_r\) chosen slots with rank-\(r\) cards: \(4!/(4-m_r)!\) ordered choices.
  2. Fill the remaining \(52 - k\) slots with the remaining \(52 - k\) cards in any order: \((52-k)!\) ways.

Dividing by \(52!\):

$$ P\!\left(\bigcap_{i \in I} A_i\right) = \frac{(52-k)!}{52!} \prod_{r=0}^{12} \frac{4!}{(4-m_r)!}. $$

This depends on \(I\) only through its multiplicity profile, and the number of size-\(k\) subsets with a given profile is \(\prod_r \binom{4}{m_r}\). So

$$ S_k = \frac{(52-k)!}{52!} \sum_{\substack{m_0+\dots+m_{12}=k \\ 0 \le m_r \le 4}} \prod_{r=0}^{12} \binom{4}{m_r} \frac{4!}{(4-m_r)!}. $$

That inner sum is a constrained sum over 13-tuples, and it is exactly the kind of thing generating functions were invented to compute.

The Generating Function

Pull out one rank's contribution. Define

$$ f(m) := \binom{4}{m} \frac{4!}{(4-m)!}, \qquad f(0), f(1), f(2), f(3), f(4) = 1,\; 16,\; 72,\; 96,\; 24. $$

Package these into a polynomial using multiplicity as the exponent of a bookkeeping variable \(x\):

$$ g(x) := \sum_{m=0}^{4} f(m)\, x^m = 1 + 16x + 72x^2 + 96x^3 + 24x^4. $$

Now expand \(g(x)^{13}\). Every term in the expansion picks one monomial \(f(m_r) x^{m_r}\) from each of the 13 factors and multiplies them together, giving \(\left(\prod_r f(m_r)\right) x^{m_0 + \dots + m_{12}}\). So the coefficient of \(x^k\) in \(g(x)^{13}\) is exactly the constrained sum we wanted:

$$ [x^k]\, g(x)^{13} = \sum_{\substack{m_0+\dots+m_{12}=k \\ 0 \le m_r \le 4}} \prod_{r=0}^{12} \binom{4}{m_r} \frac{4!}{(4-m_r)!}. $$

The upper bound \(m_r \le 4\) is enforced automatically because \(g\) has no terms above degree 4. What was a constrained nested sum over 13-tuples is now a single coefficient extraction from a fixed degree-52 polynomial.

Putting the pieces together:

$$ P(\text{win}) = \sum_{k=0}^{52} (-1)^k \frac{(52-k)!}{52!}\, [x^k]\, g(x)^{13}. $$

Sanity Check: the Single-Suit Deck

Strip the deck down to a single suit: 13 distinct cards, one per rank. Now \(f(m) = 1\) for \(m \in \{0, 1\}\), so \(g(x) = 1 + x\) and \([x^k](1 + x)^{13} = \binom{13}{k}\). The formula collapses to

$$ P(\text{win}) = \sum_{k=0}^{13} (-1)^k \binom{13}{k} \frac{(13-k)!}{13!} = \sum_{k=0}^{13} \frac{(-1)^k}{k!} \to \frac{1}{e}. $$

That is the classical derangement probability. Enumeration Solitaire with a full 52-card deck is the natural four-copy generalization: same inclusion-exclusion skeleton, but with a richer per-rank contribution polynomial than \(1 + x\).

The Answer

Before turning the crank, notice what the final formula is really gluing together. The problem mixes two different combinatorial pieces:

  • \(g(x)^{13}\) enumerates how \(k\) matches can be distributed across the 13 ranks (with each rank contributing between 0 and 4 matches).
  • \((52 - k)! / 52!\) is the probability that any particular set of \(k\) matches all happen simultaneously.

Inclusion-exclusion is what stitches them together with the alternating sign. The generating function collapses the first piece from a nested 13-tuple sum into "expand a polynomial and read off a coefficient." The second piece is a scalar weight applied per degree. So unlike the Polya change for a dollar setup, there is no single generating function whose coefficient is the answer. There is a polynomial to expand and then a 53-term weighted sum to take. That sum is the last unavoidable step.

Concretely: expand \(g(x)^{13}\) into its 53 coefficients \(c_0, c_1, \dots, c_{52}\), then compute

$$ P(\text{win}) = \sum_{k=0}^{52} (-1)^k \frac{(52-k)!}{52!}\, c_k. $$

Wolfram Alpha is not much help here. It happily expands the polynomial if you ask for

expand (1 + 16 x + 72 x^2 + 96 x^3 + 24 x^4)^13

and gives you all 53 coefficients, but its free interface will not evaluate the surrounding sum: SeriesCoefficient and Coefficient are Wolfram Language functions, and the natural-language front-end refuses to run them inside a sum[...].

Five lines of Python with sympy does the whole thing:

from math import factorial
from sympy import symbols, expand, Poly

x = symbols('x')
g = 1 + 16*x + 72*x**2 + 96*x**3 + 24*x**4
coeffs = Poly(expand(g**13), x).all_coeffs()[::-1]  # c_0, ..., c_52
p = sum((-1)**k * factorial(52 - k) * c for k, c in enumerate(coeffs)) / factorial(52)
print(float(p))  # 0.016232838...

Either route gives an exact rational with a lot of digits; the decimal value is

$$ P(\text{win}) \approx 0.01623, $$

or about 1.623%.

Two sanity checks. The Poisson heuristic \(e^{-4} \approx 0.0183\) overshoots slightly (the events are positively correlated in aggregate: once one card of a rank has appeared without matching, its remaining copies are slightly more likely to hit a later slot for that rank, which reduces the survival probability). A Monte Carlo run of 2,000,000 simulated deals gives an empirical win rate of about \(0.01633\), matching to three significant figures.

Successive games are independent, so the number of games until the first win is geometric with success probability \(\approx 0.01623\). On average you play about \(1 / 0.01623 \approx 61.6\) games before winning once. Bring a book.

What It Teaches

Enumeration Solitaire is a good showcase for two techniques that pair well:

  • Inclusion-exclusion for "probability that none of these correlated events occur." The moment you write down \(P(\bigcap A_i^c) = \sum (-1)^k S_k\), the problem becomes: what is \(S_k\)?
  • Generating functions for constrained sums over tuples. Whenever a sum runs over \((m_0, \dots, m_{n-1})\) with \(\sum m_r = k\) and per-coordinate bounds, packaging the per-coordinate contribution into a polynomial and reading off a coefficient of a product usually beats brute enumeration.

The whole solution fits on one page, and the classical derangement problem falls out as the degenerate case. Not bad for a game whose rules fit in a tweet.

References

November 2026

Baby Mac Crypto

March 2022

How to Read Ulysses

July 2020

Applied Gitflow

September 2019

Mocking AWS in Unit Tests

November 2017

A Hard(y) Math Problem