A Franciscan monk named Luca Paccioli wrote a book Summa de arithmetic, geometrica et proportionalitià in 1494, where is posed the now famous problem of points.

A and B are playing a fair game of balla. They agree to continue until one has won six rounds. The game actually stops when A has won five and B three. How should the stakes be divided.

This problem gave birth to modern probability theory. Read The Unfinished Game by Keith Devlin for a fascinating account of how Pascal and Fermat struggled to find a solution to this problem. In the book Devlin says

Interestingly, as far as we know, neither Pascal, Fermat, nor anyone else sought to resolve the issue empirically. If you were actually to play out the completion of the game many times— that is, imagine the game had been halted after eight tosses, with A ahead five to three, and then toss actual coins to complete the game—you would find that A wins roughly 7/8 of the time.

(I have reworded the above quote to match with the problem statement.)

Sounds like Devlin is asking us to verify the claim by Monte Carlo simulation. Lets do that.

import Control.Applicative ( (<$>) ) import Control.Monad.MC import Text.Printf (printf)

Lets start by defining the game. We have two players

data Player = A | B deriving (Eq, Show) players = [A,B]

Define a record to store the current state of the game.

data Game = Game { pointsA :: Int , pointsB :: Int , pointsTotal :: Int } deriving Show game :: Game game = Game 5 3 6

At any round, a player can win with equal probability. So,

play :: MC Player play = sample 2 players

Lets create a sequence of tosses. Notice that the game must end after (pointsTotal - pointsA) + pointsTotal - pointsB) - 1 rounds.

tosses :: MC [Player] tosses = sequence . replicate count $ play where count = (total - a) + (total - b) - 1 total = pointsTotal game a = pointsA game b = pointsB game

We need to keep track of the current score of each player. I use a tuple to store the score, and use the following function to update the score.

updateScore :: (Int, Int) -> Player -> (Int, Int) updateScore (a,b) A = (a+1,b) updateScore (a,b) B = (a, b+1)

Next I define a function that keeps track of the score as the game progresses.

scores :: MC [(Int,Int)] scores = scanl updateScore (a,b) <$> tosses where a = pointsA game b = pointsB game

Given a sequence of scores, I can determine who won the game.

winner :: [(Int, Int)] -> Player winner [] = error "Calculating winner of an empty sequence" winner ((a,b) : xs) | a >= total = A | b >= total = B | otherwise = winner xs where total = pointsTotal game

Now, we need to compute the probability that A wins the game. Since, I am using a Monte Carlo simulation, I need a random variable to keep track of the number of times player A has won. This can be done by simply assigning a reward of 1 when player A wins and a reward of 0 when player 2 wins. That way, the expected value of the reward will be equal to the probability that A wins.

reward :: Player -> Double reward A = 1 reward B = 0

And finally, the Monte Carlo simulation.

main :: IO () main = let n = 1000000 seed = 101 stats = repeatMC n (reward <$> winner <$> scores) `evalMC` mt19937 seed in do printf "Probability that A wins : %f

" (sampleMean stats) printf "99%% Confidence interval : (%f, %f)

" `uncurry` (sampleCI 0.99 stats)

This gives

Probability that A wins : 0.8749600000000389 99% Confidence interval : (0.8741080072900288, 0.8758119927100491)

The actual solution of the problem is 7/8=0.875 . So, the problem would have been resolved much earlier had the Renascence mathematician had access to Monte Carlo simulations 🙂