Keith Devlin's book The Unfinished Game describes how Fermat and Pascal discovered the rules of probability that guide gambling in games. The question they confronted was: what if a gambling game is interrupted, but one player is in the lead by a certain score. How much of the pot should the leader get?
My friends and I faced a similar question when a game of Risk ran on too long (as they often do). Player A was poised to make a sweeping attack on player D, whose territories were arranged in such a way that A could attack from one territory to the next without ever having to branch off. We wrote down the number of A's massed armies, 72, and the number of armies in D's successive territories: 22, 8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1. What is the probability that A can capture all these territories?
Let's explain some Risk rules and terminology:
A battle is when armies from one territory attack an enemy neighboring territory. A roll of dice determines which armies perish:
A campaign consists of a sequence of battles and occupations. In this notebook we will consider only a chain campaign, in which attackers invade successive enemy territories in order, always moving all their remaining armies into each captured territory, never branching, and never changing strategy based on the outcome of a battle. The attackers win the campaign if they occupy all the territories. The attackers lose if there are any remaining defenders and only one remaining attacker, who by rule cannot attack.
With that out of the way, we're ready for some Python implementation.
from typing import List, Iterable, Tuple
from collections import Counter
from functools import lru_cache
from dataclasses import dataclass
import random
import itertools
import matplotlib.pyplot as plt
I will represent the state of a campaign with the class State
, and the state of the unfinished game as start
:
@dataclass(frozen=True)
class State:
A: int # Number of attackers
D: int # Number of defenders in first territory
rest: Tuple[int] = () # Tuple of numbers of defenders in subsequent territories
start = State(A=72, D=22, rest=(8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1))
The function update
will update a state to reflect the Deaths
that happened in a battle. We declare game_over
when there are no defenders or only a single attacker (who can't leave their territory) remaining.
@dataclass(frozen=True)
class Deaths:
"The number of attackers and defenders who die in a battle."
A: int
D: int
def update(state: State, dead: Deaths) -> State:
"""Update the `state` of a campaign to reflect the`dead` in a battle."""
a = state.A - dead.A # Attackers remaining
d = state.D - dead.D # First territory defenders remaining
r = state.rest # Other territories defenders remaining
return (State(a, d, r) if d # Defenders still in first territory
else State(a - 1, r[0], r[1:]) if r # First territory captured
else State(a - 1, 0)) # All territories captured
def game_over(state) -> bool:
"""Is the game over?"""
return state.D == 0 or state.A <= 1
We'll represent a roll of the dice with a list of integers; for example the attacker might roll [2, 6, 1]
with three dice. The function random_roll(n)
gives a random outcome for rolling n
dice. The function battle_deaths
, when given the two lists of dice rolled by attacker and defender, returns the number of attackers and defenders who perish in the battle.
die = (1, 2, 3, 4, 5, 6)
Dice = List[int] # a list of die rolls, like [2, 6, 1]
def random_roll(n) -> Dice:
"""Roll n dice randomly."""
return [random.choice(die) for _ in range(n)]
def battle_deaths(A_dice: Dice, D_dice: Dice) -> Deaths:
"""How many (attacker, defender) armies perish as the result of these dice?"""
dead = Counter('D' if a > d else 'A'
for a, d in zip(sorted(A_dice, reverse=True),
sorted(D_dice, reverse=True)))
return Deaths(dead['A'], dead['D'])
# Example battle
A_dice = random_roll(3)
D_dice = random_roll(2)
A_dice, D_dice, battle_deaths(A_dice, D_dice)
([6, 1, 6], [6, 4], Deaths(A=1, D=1))
A simulation makes random choices at every choice point (here, every dice roll), and reports on the outcome. The function simulate_campaign
rolls random dice for a battle until the game is over:
def simulate_campaign(state) -> State:
"""Simulate a campaign with random rolls, returning the final state."""
while not game_over(state):
dead = battle_deaths(random_roll(min(3, state.A - 1)),
random_roll(min(2, state.D)))
state = update(state, dead)
return state
Let's see who wins:
simulate_campaign(start)
State(A=5, D=0, rest=())
That final state says that the attackers won—there are 5 attackers and no defenders left.
But what if we want to know the probability that the attackers win? Monte Carlo simulation answers that question by repeating a simulation many times and summarizing all the final outcomes. The summary is in the form of a probability distribution, a mapping of {outcome_state: probability}
pairs, which I have implemented as ProbDist
.
Note: I have ProbDist
as a subclass of Counter
, with the restriction that the values are normalized to sum to 1. I realize that the name Counter
suggests integer counts, but that's not a requirement, and Counter
has a nicer API than dict
.
class ProbDist(Counter):
"A Probability Distribution."
def __init__(self, *args):
"Normalize total to 1."
Counter.__init__(self, *args)
total = sum(self.values())
for x in self:
self[x] /= total
ProbDist(random_roll(1000)) # Probability distribution over 1,000 die rolls
ProbDist({4: 0.176, 3: 0.168, 2: 0.163, 1: 0.175, 5: 0.16, 6: 0.158})
The higher-order function monte_carlo
works with any random simulation function, calling it k
times, each time passing in a start
state, and collecting the k
outcome states in a ProbDist
:
def monte_carlo(simulation, start, k=1000) -> ProbDist:
"Call `simulation(start)` repeatedly (`k` times) and return a ProbDist of outcomes."
return ProbDist(simulation(start) for _ in range(k))
Here we simulate the campaign 10 times and see the summary of outcomes:
monte_carlo(simulate_campaign, start, 10)
ProbDist({State(A=16, D=0, rest=()): 0.1, State(A=1, D=5, rest=(1,)): 0.1, State(A=11, D=0, rest=()): 0.2, State(A=1, D=2, rest=(5, 1)): 0.1, State(A=17, D=0, rest=()): 0.1, State(A=6, D=0, rest=()): 0.2, State(A=1, D=0, rest=()): 0.1, State(A=28, D=0, rest=()): 0.1})
And here we run 1,000 simulations and report how often the attackers win:
def attacker_win_probability(dist: ProbDist) -> float:
"""The probability that the attackers win the campaign."""
return sum(dist[s] for s in dist if not s.D)
attacker_win_probability(monte_carlo(simulate_campaign, start, 1000))
0.8190000000000004
The Monte Carlo simulation says the attackers win about 82% of the time. How accurate is that result? The standard deviation of the expected value of a binomial variable is $\sqrt{p(1-p)/n}$, where $p$ is the true probability of one of the two outcomes and $n$ is the number of samples. In our case $\sqrt{0.8(1-0.8)/1000}$ gives a standard deviation of about 1%. So we can be pretty confident that the true percentage is within 3 standard deviations; that is, between 79% and 85%.
We could get better accuracy at the cost of increased computing time. To get the standard deviation down by a factor of 10 from 1% to 0.1% would require 100 times more computation (because of the square root in the formula).
An alternative to the Monte Carlo approach is to explicitly calculate the exact probability distribution of the final state of the campaign. The differences between the two approaches are:
The Monte Carlo approach deals with a single current state, using random dice rolls to decide how that state changes. At the end of a simulated campaign we get a single final state, and we repeat the simulation to get an estimated probability distribution over final states.
The exact probability calculation approach deals with a probability distribution over states, right from the start. At each dice roll, every possible outcome is considered, and the probability distribution is updated to reflect all the possible outcomes and their probabilities. At the end of the campaign we have an exact probability distribution over all possible final states (well, exact at least to the limits of floating point precision; if you need more precision, use fractions.Fraction
).
In deciding whether to use the exact calculation approach, there are three considerations:
Let's examine the three considerations:
Worth the effort? If I only wanted to know whether the attackers chance of winning is above or below 50%, then the Monte Carlo approach is the quickest way to answer the question. The code will be simple and straightforward. Monte Carlo code deals with the single current state of the simulation, for example:
while not game_over(state):
But in the exact calculation we need to consider every possible states (referenced in a loop):
while any(not game_over(state) for state in probdist):
Putting all these loops into the code makes it more complex. (Less added complexity in Python, which has comprehensions, than in other languages where loops must be statements, not expressions.) After the shared basics of representing states and battle outcomes, we need just 10 non-comment lines of code to implement the Monte Carlo method (in simulate_campaign
, monte_carlo
, and random_roll
above), but twice as much code (22 lines) to implement exact calculation (in battle_deaths_probdist
, all_rolls
, campaign_probdist
, and battle_update_probdist
below).
Possible? Imagine a Risk rule change where if the attackers roll three 6s and the defenders roll two 6s, then both sides gain an army. For the Monte Carlo simulation it would be trivial to add a single if
statement to handle this rule, and the expected run time of the program would hardly change. But with exact calculation everything has changed, because we now have an infinite game: no matter how many moves ahead we look, there will always be some possible states that have not terminated. (To deal with these we would have to make some compromises, such as stopping the calculation when there are still some nonterminal states, as long as they have sufficiently low total probability. Sometimes a mathematical formula can determine the value of an infinite game.)
Feasible? Imagine a rule change where before each battle, every army independently has the option to take one of ten possible actions (e.g. to move to a safe neighboring territory). Then with just 80 armies, the very first move has $10^{80}$ possible outcomes, meaning that the probability distribution requires as many states as there are atoms in the universe. (To deal with this we would either stick with the Monte Carlo simulation, or a variant of it such as particle filtering in which we maintain a sample of several possible states–more than a single state, but less than the complete state distribution, or we would find a way to abstract over the possible moves.)
The Risk campaign problem as it stands leads to a very efficient exact probability calculation (possible, feasible, and, IMHO, worth it). If there are $n$ total armies to start, then there are fewer than $n^2$ possible states, and the game can last no more than $n$ moves. With $n$ in the range of a few hundred, computation takes less than a second; much faster and more accurate than doing 100,000 simulations.
The function battle_deaths
was defined above to return the specific death counts for a specific dice roll.
Now we'll define battle_deaths_probdist
to give a probability distribution over all possible battle outcomes, corresponding to all possible dice rolls. The input to battle_deaths_probdist
is the number of attacking and defending armies for just this one battle (not the total number of attacking and defending armies in the current state). Thus, the attackers will always be 3, 2, or 1, and the defenders will always be 2 or 1. I define it this way so that I can cache the results and reuse them in subsequent battles, for efficiency.
@lru_cache()
def battle_deaths_probdist(battlers: State) -> ProbDist:
"""A probability distribution of deaths in a single battle.
Requires 1 <= battlers.A <= 3 and 1 <= battlers.D <= 2."""
return ProbDist(battle_deaths(A_dice, D_dice)
for A_dice in all_rolls(battlers.A)
for D_dice in all_rolls(battlers.D))
def all_rolls(n) -> Iterable[Dice]:
"""All possible rolls of `n` dice."""
return tuple(itertools.product(die, repeat=n))
battle_deaths_probdist(State(3, 2)) # Three attacker dice against two defender dice
ProbDist({Deaths(A=2, D=0): 0.2925668724279835, Deaths(A=1, D=1): 0.3357767489711934, Deaths(A=0, D=2): 0.37165637860082307})
Our old function campaign
is mimiced by campaign_probdist
, except that the former updates a State
whereas the later updates a ProbDist
. We start with certainty: there is a 100% chance that we are initially in the start
state. But the fog of war means that uncertainty soon arises: we don't know how the dice will land, so we don't know the outcome of the very first battle. Subsequent battles add more uncertainty.
campaign_probdist
calls battle_update_probdist
to update the probability distribution to account for one battle, in all the possible ways it can turn out. The key line in battle_update_probdist
is
outcomes[update(state, dead)] += probdist[state] * dead_probdist[dead]
which tells us to consider the outcome in which dead
is the number of attackers and defenders that die in a battle, and to update state
with those death counts, and increment the probability of that updated outcome by the product of the probability of state
and the probability of dead
given state
.
def campaign_probdist(start: State) -> ProbDist:
"""Probability distribution for all outcomes of a campaign."""
probdist = ProbDist({start: 1.0})
while any(not game_over(state) for state in probdist):
probdist = battle_update_probdist(probdist)
return probdist
def battle_update_probdist(probdist) -> ProbDist:
"""For every possible campaign state in the `probdist`, consider the outcomes of a battle.
Combine these all into one updated `outcomes` ProbDist."""
outcomes = ProbDist()
for state in probdist:
if game_over(state): # `state` carries through unchanged to `outcomes`
outcomes[state] += probdist[state]
else: # Replace `state` with all possible outcomes from a battle
dead_probdist = battle_deaths_probdist(Deaths(min(3, state.A - 1), min(2, state.D)))
for dead in dead_probdist:
outcomes[update(state, dead)] += probdist[state] * dead_probdist[dead]
return outcomes
What is the exact probability of winning the unfinished game?
attacker_win_probability(campaign_probdist(start))
0.8105485936352178
The attackers defeat all the defenders 81.05% of the time.
What are the 20 most common outcomes?
campaign_probdist(start).most_common(20)
[(State(A=12, D=0, rest=()), 0.03824220182706657), (State(A=11, D=0, rest=()), 0.0380992150215239), (State(A=13, D=0, rest=()), 0.038032667992797725), (State(A=10, D=0, rest=()), 0.037618411457469664), (State(A=14, D=0, rest=()), 0.03746512036620402), (State(A=9, D=0, rest=()), 0.036822601595588304), (State(A=15, D=0, rest=()), 0.036544033638847624), (State(A=8, D=0, rest=()), 0.035741340182918115), (State(A=16, D=0, rest=()), 0.035284184613703425), (State(A=7, D=0, rest=()), 0.03440940023374102), (State(A=17, D=0, rest=()), 0.03371050842173721), (State(A=6, D=0, rest=()), 0.03286517629527088), (State(A=18, D=0, rest=()), 0.031857448871758426), (State(A=5, D=0, rest=()), 0.031149102741286083), (State(A=19, D=0, rest=()), 0.029767807751056283), (State(A=4, D=0, rest=()), 0.029302161098148583), (State(A=20, D=0, rest=()), 0.027491133463900305), (State(A=3, D=0, rest=()), 0.0273645356924103), (State(A=1, D=5, rest=(1,)), 0.026621109937678585), (State(A=1, D=1, rest=()), 0.025661022694069335)]
The most probable outcome is that the attackers win with 12 remaining armies. The top 10 outcomes have anywhere from 7 to 16 attackers remaining. You have to go down to the 19th most common outcome to find one where the defenders win.
Let's try to visualize all the possible outcomes. I'll define the score of a campaign as the number of attacker armies in the final territory minus the total number of defenders, minus the number of territories that the defenders hold. This score will be positive when the attackers win and negative when they lose. A score cannot be zero. The function show
plots score probabilities, with the attacker's wins in green and the defender's wins in blue:
def show(probdist, epsilon=0.0002):
"""Plot and annotate a probability distribution over states."""
states = [s for s in probdist if probdist[s] > epsilon] # Ignore low-probability states
X = [score(s) for s in states]
Y = [probdist[s] for s in states]
μ = sum(score(s) * probdist[s] for s in probdist)
p = attacker_win_probability(probdist)
plt.figure(figsize=(10, 5))
plt.title(f'Attacker wins {p:.2%}. Average score: {μ:.2f}.')
plt.xlabel('Score (positive when attacker wins)')
plt.ylabel('Probability of score')
plt.bar(X, Y, width=0.7, color=['g' if x > 0 else 'b' for x in X])
def score(state):
return (state.A if not state.D else
state.A - (state.D + sum(state.rest)) - (len(state.rest) + 1))
show(campaign_probdist(start))
Interesting! To the right of 0, we see a nice bell-shaped curve. But there is a decidely non-normal pattern on the left side. There are gaps—scores for which the probability is zero—and there are spikes that rise above the normal curve.
What's causing the gaps? We know that a campaign score can never be 0. And in this campaign, no score can be -2, because a -2 can only be achieved with the final state State(1, 2)
: one attacker and two defenders in a single territory; this campaign has only one defender in the final territory, so State(1, 2)
could never occur. Looking at the plot, we see gaps surrounding groups of bars, where the group sizes, reading left to right, are 1, 1, 3, 1, 2, 3, 5, and 1. These are exactly the number of defenders in the final 8 territories, so the gaps are indicating the impossible scores.
What's causing the spikes? To some extent they represent probability mass that has shifted over a place from the impossible states to the neighboring possible states. But the spikes are not quite tall enough to fill in the gaps. I admit I don't understand the exact shape of each spike–do you?
To help your intuition, here's a campaign where 100 attackers take on 100 defenders, 50 in the first territory and 5 in each of 10 more territories:
show(campaign_probdist(State(100, 50, 10 * (5,))))
In the plot below all the defenders are in a single territory, and there is no visible gap–the bell shape is restored. However, there is something going on with odd-versus-even scores. I believe that part of what is happening is that all but the very last few battles will be 3-versus-2 battles, and thus will result in 2 casualties, preserving the parity of the total number of armies. The parity is broken only in the final battles. But I don't think that's the full story—what do you think?
show(campaign_probdist(State(146, 100)))
The odd-even disparity can be lessened by adding a few territories where there will not be 3-versus-2 battles. The result is a very smooth curve, although not quite normal—it is clearly a bit fatter on the left than the right (the mean is 60 and the probability of 20 is greater than the probability of 100.
show(campaign_probdist(State(146, 93, (1, 1, 1, 1))))
I want to analyze a simpler situation: With A attackers and D defenders in a single territory, what is the probability the attackers win?
I will make a chart with the number of defenders varying from 1 to 60, and with the number of attackers separated into eight cases (depicted as eight lines), where in each case there are Δ more attackers than defenders:
def chart(defenders=range(1, 61), deltas=(9, 5, 2, 1, 0, -1, -2, -5)):
"""Plot win probability for various numbers of defenders and Δ more attackers."""
plt.figure(figsize=(10, 5)); plt.grid()
plt.title('Each line: attackers with Δ more armies than defenders')
plt.xlabel('Number of Defenders'); plt.ylabel('Win Probability for Attackers')
plt.yticks([i / 10 for i in range(11)]); plt.xticks(range(0, 61, 5))
for delta in deltas:
winprobs = [attacker_win_probability(campaign_probdist(State(max(0, d + delta), d,)))
for d in defenders]
plt.plot(defenders, winprobs, 'o-', label=f'Δ = {delta:2}')
plt.legend(shadow=True)
chart()
Note that the purple line (fourth from bottom), where the number of attackers exactly equals the number of defenders, gives a low win probability for a small attacking force, but reaches 50% for 12-on-12, and 73% for 60-on-60. The red line, where the attackers have one more army than the defenders, dips from one to two defenders but is over 50% for a 6-on-5 attack. Similarly, the green line, where the attackers have a surplus of two armies, dips sharply from 75% to 66% as the number of defenders goes from 1 to 2, dips slightly more for 3 and 4 defenders, and then starts to rise. So overall, an attacker does not need a big advantage in armies as long as there are many armies on both sides. Even when the attacker is at a disadvantage in numbers (as in the bottom grey line where the attacker has five fewer armies), the attacker still has better than 50/50 odds with 45 attackers or more.
Here are regression unit tests that also serve as examples of usage.
# Tests for States
assert start.A == 72
assert start.D == 22
assert start.rest == (8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1)
assert sum(start.rest) + start.D == 60
assert len(start.rest) == 13
# Tests for update
small = State(10, 2, (3, 4)) # Small example where 10 attackers battle 2 + 3 + 4 defenders
assert update(small, Deaths(0, 2)) == State(9, 3, (4,))
assert update(small, Deaths(1, 1)) == State(9, 1, (3, 4))
assert update(small, Deaths(2, 0)) == State(8, 2, (3, 4)) # continue to battle
assert update(State(10, 2, ()), Deaths(0, 2)) == State(9, 0, ())
# The 4 sample dice rolls from www.ultraboardgames.com/risk/game-rules.php
assert battle_deaths([1, 1, 6], [3,]) == Deaths(0, 1) # Defender loses one
assert battle_deaths([2, 6, 1], [2, 3]) == Deaths(1, 1) # Both lose one
assert battle_deaths([3, 3], [3, 4]) == Deaths(2, 0) # Attacker loses two
assert battle_deaths([6,], [5, 4]) == Deaths(0, 1) # Defender loses one
# All six possible battle dice combinations. The answers match those given at
# a Risk data analysis blog: http://datagenetics.com/blog/november22011/
assert battle_deaths_probdist(State(1, 1)) == ProbDist(
{Deaths(1, 0): 21, Deaths(0, 1): 15})
assert battle_deaths_probdist(State(2, 1)) == ProbDist(
{Deaths(1, 0): 91, Deaths(0, 1): 125})
assert battle_deaths_probdist(State(3, 1)) == ProbDist(
{Deaths(1, 0): 441, Deaths(0, 1): 855})
assert battle_deaths_probdist(State(1, 2)) == ProbDist(
{Deaths(1, 0): 161, Deaths(0, 1): 55})
assert battle_deaths_probdist(State(2, 2)) == ProbDist(
{Deaths(2, 0): 581, Deaths(1, 1): 420, Deaths(0, 2): 295})
assert battle_deaths_probdist(State(3, 2)) == ProbDist(
{Deaths(2, 0): 2275, Deaths(1, 1): 2611, Deaths(0, 2): 2890})
assert set(all_rolls(2)) == {
(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6),
(2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6),
(3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6),
(4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6),
(5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6),
(6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6)}
assert attacker_win_probability(campaign_probdist(State(2, 1))) == 15/36
assert campaign_probdist(State(2, 1)) == ProbDist({State(1, 1): 21/36, State(1, 0): 15/36})
assert campaign_probdist(State(4, 2, (1,))) == ProbDist(
{State(1, 2, (1,)): 0.21807067805974698,
State(1, 1, (1,)): 0.12597532213712342,
State(1, 1): 0.294669783648961,
State(1, 0): 0.14620529335276633,
State(2, 0): 0.21507892280140226})
assert battle_update_probdist(ProbDist({start: 1})) == ProbDist(
{State(70, 22, (8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1)): 0.2925668724279835,
State(71, 21, (8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1)): 0.3357767489711934,
State(72, 20, (8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1)): 0.37165637860082307})
True
True