# Problem of Points using Monte Carlo

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 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) 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\n" (sampleMean stats)
printf "99%% Confidence interval : (%f, %f)\n"
`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 :-)