The 538 Riddler Classic for 28 August 2020 asks about the chance of winning the children's card game War in a sweep: a game where player A wins 26 turns in a row against player B. (In War, players are dealt 26 cards each and on each turn they flip the top card in their respective hands; higher rank card wins. Other rules are not relevant to this problem.)
We'll analyze this problem and come up with a program to solve it. First let's get the imports out of the way:
import random
import collections
from itertools import permutations, combinations
from statistics import mean
Here are the four approaches I considered:
There is no strategy in War, so on every turn, both players have an equal chance of winning the turn. If player A wins each turn with probability 1/2, than the probability of sweeping would be $(1/2)^{26}$ or about 15 in a billion. But that's wrong because there is also the possibility of a tie. After A picks a card, there are 51 cards remaining, and 3 have the same rank as A's card, so the possible outcomes for the first turn are:
tie: 3/51
A wins: 24/51
B wins: 24/51
If every subsequent turn were the same, the probability that A sweeps would be $(24/51)^{26}$ or about 3 in a billion. But that's not quite right because the outcome of each turn depends on the cards picked in the previous turns. So simple arithmetic isn't enough.
Brute force enumeration means:
Easy-peasy; here's the code:
Probability = float # Type: Probability is a number between 0.0 and 1.0
Deck = list # Type: Deck is a list of card ranks
def make_deck(ranks=13, suits=4) -> Deck:
"""A Deck is a list of ranks, each rank repeated `suits` times."""
return [r + 1 for r in range(ranks) for _ in range(suits)]
def sweep(deck) -> bool:
"""Upon dealing this deck, does player A win every turn?"""
return all(deck[i] > deck[i + 1] for i in range(0, len(deck), 2))
def p_sweeps(decks) -> Probability:
"""The probability that A sweeps, averaged over the decks."""
return mean(sweep(deck) for deck in decks)
(Implementation details: in Python False
is equivalent to 0
and True
is equivalent to 1
, so the mean
of an iterable of bool
values is the same as the proportion (or probability) of True
values. Also, it just felt wrong to have a card with rank 0
, thus I add 1 so that ranks start at 1
.)
Here are two 8-card decks:
make_deck(8, 1)
[1, 2, 3, 4, 5, 6, 7, 8]
make_deck(4, 2)
[1, 1, 2, 2, 3, 3, 4, 4]
And here are the probabilities of sweeping a game played with those decks:
p_sweeps(permutations(make_deck(8, 1)))
0.0625
p_sweeps(permutations(make_deck(4, 2)))
0.03571428571428571
What about the real deck, with 52 cards? Unfortunately, there are 52! permutations (more than $10^{67}$), and even if we were clever about the duplicated ranks and the ordering of the 26 turns, and even if we could process a billion deals a second, it would still take millions of years to complete the brute force enumeration. And 538 Riddler wanted the answer by Monday.
It takes too long to look at all the permutations, but we can randomly sample the permutations. We call that a simulation:
def sample(deck, N) -> Deck:
"""Yield N randomly shuffled deals of deck."""
for _ in range(N):
random.shuffle(deck)
yield deck
p_sweeps(sample(make_deck(), 10_000))
0
Sampling 10,000 deals with the full 52-card deck we got zero sweeps. Approach 1 suggests the probability of a sweep is somewhere near 3 in a billion, which means we would need to sample trillions of deals to get a reliable estimate of the probability. That would require roughly a day of run time: much better than millions of years, but still not good enough.
As discussed in my How To Count Things notebook, in problems where brute force enumeration is not feasible it is often possible to consider every outcome if we use a representation that is:
For example, if there are two cards remaining then there are 52 × 51 possible concrete decks, but only two abstract cases:
[5, 5]
or [9, 9]
, or any other pair of cards.[10, 8]
or [2, 5]
or any other two distinct ranks.The representation I will use, which I call AbDeck
for abstract deck, is as follows:
AbDeck
is a tuple of counts of the number of cards of each rank, without specifying the rank.(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
(1, 1)
and (2,)
.(1, 1, 1, 1)
, (1, 1, 2)
, (1, 3)
, (2, 2)
, and (4,)
.AbDeck = tuple # Type: An abstract deck
def make_abdeck(ranks=13, suits=4) -> AbDeck:
"""Make an abstract deck."""
return AbDeck([suits] * ranks)
deck = make_abdeck(13, 4) # The initial 52-card abstract deck
deck
(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)
*Note:* I originally thought of two possible abstract representations for decks:
AbDeck
as described above, and an alternative in which:
deck[i]
gives the number of different ranks that have exactly i
cards remaining in the deck.(0, 0, 0, 0, 13)
: 13 ranks with 4-of-a-kind.(0, 2)
(two ranks with 1 card each) and (0, 0, 1)
(one rank with two cards).That representation was more compact, and ran faster; I used it in the first version of this notebook. But the very nice implementation shared by Laurent Lessard convinced me to switch to the
AbDeck
representation becausse it is simpler to code and understand.
We want to keep track of the possible outcomes of each turn. We'll do that with the help of a class called PDist
(for probability distribution), a mapping of {deck: p, ...}
where deck
is an abstract deck that is the outcome after some number of turns, each of which A has won, and p
is the probability of that outcome.
class PDist(collections.Counter):
"""A probability distribution of {abstract_deck: probability}."""
The probability that A sweeps a game of War given an abstract deck, p_sweep(deck)
, is computed as follows:
P
being a probability distribution of outcomes after 0 turns: {deck: 1.0}
.P
to be the probability distribution that results from playing a turn, P = play_turn(P)
:deck
in P
:i
and rank j
from deck
to yield a resulting deck:deck
× probability of i
× probability of j
).P
from the final turnNote that (for example) combinations(range(3), 2)
yields (1, 2)
, but not (1, 1)
nor (2, 1)
. That's just what we want: if the two selected cards are i, j = (1, 1)
then the result will be a tie, which we don't want to put in the probability distribution of decks where A sweeps. And out of the two possibilities i, j = (1, 2)
and i, j = (2, 1)
, exactly one will result in a win for A, so we only want to include one of them.
def p_sweep(deck: AbDeck) -> Probability:
"""The probability that player A sweeps a game of War."""
P = PDist({deck: 1.0}) # The initial probability distribution
for turn in range(sum(deck) // 2):
P = play_turn(P)
return sum(P.values())
def play_turn(P) -> PDist:
"""Play one turn with all possible card choices for players A and B; return
the probability distribution of outcomes where A still might sweep."""
P1 = PDist() # The probability distribution of {deck: prob} after this turn
for deck in P:
n = sum(deck) # number of cards in deck
for i, j in combinations(range(len(deck)), 2):
P1[remove(deck, i, j)] += P[deck] * deck[i] / n * deck[j] / (n - 1)
return P1
def remove(deck, *indexes) -> AbDeck:
"""Remove cards at indexes from deck."""
counts = [c - indexes.count(i) for i, c in enumerate(deck)]
return AbDeck(sorted(c for c in counts if c > 0))
What's the probability that player A will win in a sweep?
p_sweep(deck)
3.132436174322299e-09
A little over 3 in a billion.
(By the way, this computation took about 0.2 seconds, a big improvement over millions of years.)
That's the answer to my question about the probability of A sweeping, but 538 Riddler actually posed the question somewhat differently: "How many games of War would you expect to play until you had a game that lasted just 26 turns with no wars, like Duane’s friend’s granddaughter?" That is, they want to know the inverse of the probability, and they are allowing for either A or B to sweep. So that would be:
1 / (2 * p_sweep(deck))
159620171.70491105
We would expect a sweep about once every 160 million games. I have to say, I'm begining to doubt Duane’s friend’s granddaughter's claim.
Let's work through how p_sweep
and play_turn
work. We'll start with P
being the initial probability distribution, and then play four turns, each time displaying how P
is updated:
P = PDist({deck: 1.0})
P
PDist({(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 1.0})
P = play_turn(P) # Turn 1
P
PDist({(3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.47058823529411703})
That makes sense: the only result where A still has a chance to sweep is when the two players select different ranks, giving us two ranks with three-of-a-kind remaining (and 11 with four-of-a-kind remaining). The probability of A sweeping one turn is 24/51 = 0.47058823529411703.
P = play_turn(P) # Turn 2
P
PDist({(2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.001728691476590634, (2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.050708283313325254, (3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.16902761104441777})
P = play_turn(P) # Turn 3
P
PDist({(1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 3.065055809557862e-06, (1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.00040458736686163773, (2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.010114684171540949, (1, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.0017981660749406122, (2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.00020229368343081884, (2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4): 0.04855048402339655, (3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4): 0.043155985798574735})
P = play_turn(P) # Turn 4
P
PDist({(4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 1.4807032896414793e-09, (1, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 5.212075579538007e-07, (1, 1, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 3.648452905676605e-05, (2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 5.863585026980257e-07, (2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 1.563622673861402e-05, (1, 1, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4): 2.3454340107921024e-06, (1, 2, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4): 0.00018763472086336856, (1, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4): 0.0016887124877703135, (2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4): 4.397688770235195e-05, (2, 2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4): 0.0023923426910079496, (2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4): 0.014635508227342713, (3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4): 3.127245347722804e-05, (1, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4): 0.002001437022542595, (2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4): 0.02101508873669728, (3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4): 0.007005029578899085})
We'll leave it as an exercise for the reader to verify that these turns are correct.
One way to gain confidence: The answer should be somewhere close to the arithmetic calculation of $(24/51)^{26}$.
p_sweep(deck) / ((24/51) ** 26)
1.0167508045532485
We see that it is indeed close, differing by less than 2%.
Another way to gain confidence: We can verify that there is no difference (∆ = 0) between the brute force p_sweeps
and the abstract incremental p_sweep
for small decks. This gives us some confidence that both are correct (or maybe both have similar errors). p_sweeps
takes less than 50 milliseconds to solve decks of 8 cards or less, but takes about 4 seconds to do the 10! (almost 4 million) permutations of ten card decks.
for r, s in [(2, 1),
(4, 1), (2, 2),
(6, 1), (3, 2), (2, 3),
(8, 1), (4, 2), (2, 4),
(5, 2), (2, 5)]:
p1 = p_sweeps(permutations(make_deck(r, s)))
p2 = p_sweep(make_abdeck(r, s))
print(f'deck({r}, {s}): p = {p1:.4f}; ∆ = {abs(p1 - p2):.9f}')
deck(2, 1): p = 0.5000; ∆ = 0.000000000 deck(4, 1): p = 0.2500; ∆ = 0.000000000 deck(2, 2): p = 0.1667; ∆ = 0.000000000 deck(6, 1): p = 0.1250; ∆ = 0.000000000 deck(3, 2): p = 0.0667; ∆ = 0.000000000 deck(2, 3): p = 0.0500; ∆ = 0.000000000 deck(8, 1): p = 0.0625; ∆ = 0.000000000 deck(4, 2): p = 0.0357; ∆ = 0.000000000 deck(2, 4): p = 0.0143; ∆ = 0.000000000 deck(5, 2): p = 0.0180; ∆ = 0.000000000 deck(2, 5): p = 0.0040; ∆ = 0.000000000