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 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\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 🙂

Advertisements

Monty Hall with a billion doors

In order to provide intuition behind the solution of the Monty Hall problem, Antonio Cangiano says:

If there were a billion doors, you picked one, and then Monty proceeded to open up all the remaining doors but one, we’d have a situation where it would be extremely unlikely that you picked the right door at the beginning, while it would be extremely likely that the remaining door was the one that was concealing the car.

Now, what happens when there are more than three doors. In this post, I will modify the solution of my last post to work for any number of doors. It requires a little change to the program.

 import Data.List (delete, (\\) )
 import Text.Printf (printf)
 import Control.Monad
 import Control.Monad.MC

This time, instead of defining Door by listing all alternatives, I define it as an instance of Bounded class.

 data Door   = Door Int deriving (Show, Eq)
 numDoors    = 3

 instance Bounded Door where
   minBound = Door 1
   maxBound = Door numDoors

 instance Enum Door where
   toEnum      n     = Door n
   fromEnum (Door n) =    n

 doors = [minBound .. maxBound] :: [Door]

Object are the same as before.

 data Object = Car | Goat deriving (Show, Eq)
 type Universe = [Object]

There is one car; all other doors have goats. This can still be generated using the shuffle function.

 nature :: MC Universe
 nature = shuffle numDoors (Car: repeat Goat)

The contestant open one door at random.

 contestant :: MC Door
 contestant = sample numDoors doors

The options available to the host remain the same as before.

 options :: Universe -> Door -> [Door]
 options universe door = filter (\d -> Goat == open universe d) remainder
   where remainder = delete door doors

We need a generic function to see what is behind a closed door.

 open :: Universe -> Door -> Object
 open universe (Door n) = universe !! (n-1)

The host opens all but one door. We can do this by replacing the sample function with sampleSubset function.

 host :: [Door]-> MC [Door]
 host o = sampleSubset (numDoors-2) l o where l = length o

The strategies are now whether to stick or switch based on the door that the contestant chose, and the doors that the host opened.

 strategyStick :: Door -> [Door] -> Door
 strategyStick d1 d2 = d1

 strategySwitch :: Door -> [Door] -> Door
 strategySwitch d1 d2 = head (doors \\ (d1:d2))

The reward function is same as before.

 reward :: Universe -> Door -> Double
 reward u d = if Car == open u d then 1.0 else 0.0

And so is the generation of the ranom experiment.

 randomVariable :: (Door -> [Door] -> Door) -> MC Double
 randomVariable strategy = do
   u <- nature
   c <- contestant
   h <- host (options u c)
   let r = reward u (strategy c h)
   return r

And we see the result.

 main = let
   n     = 100000
   seed  = 42
   stats = repeatMC n (randomVariable strategyStick) `evalMC` mt19937 seed
   in do
      putStrLn $ printf "\nnumDoors                : %d" numDoors
      putStrLn $ printf "Mean                    : %f" (sampleMean stats)
      putStrLn $ printf "99%% Confidence interval : (%f, %f)"
                `uncurry` (sampleCI 0.99 stats)

Now lets see the results. For

numDoors                : 3
Mean                    : 0.3334800000000033
99% Confidence interval : (0.3296397390134338, 0.33732026098657275)

numDoors                : 30
Mean                    : 0.03466999999999957
99% Confidence interval : (0.03317983597888913, 0.03616016402111001)

numDoors                : 300
Mean                    : 0.0038399999999999784
99% Confidence interval : (0.0033362101507491866, 0.00434378984925077)

Well, I think that you get the idea. I am not going to run this thing for a billion doors.

Monty Hall problem using Monte Carlo simulations

n0ne showed how to translate Monte Carlo simulations of the Monty Hall Problem from Ruby and Python to Haskell. Since I have been trying to understand Patrick Perry’s monte-carlo monad, I thought that this is a good problem to start with. Monty Hall problems is one of those probability problems where the result is opposite to what intuition suggests.

Lets start by importing some Haskell packages that I will use later on. The most important for our purpose is Control.Monad.MC, which provides a nice wrapper to do Monte Carlo simulations.

 import Data.List (delete, (\\) )
 import Text.Printf (printf)
 import Control.Monad
 import Control.Monad.MC

So, the Monty Hall problem goes as follows. A contestant in game show faces three doors. There is a car behind one of them and goats behind the other two. He picks a door, and then the game show host, who knows what is behind the doors, opens one of the other doors to reveal a goat. The contestant then has the option to stick to his original choice, or switch to the unopened door. What is better?

One thing that I like about Haskell is its expressiveness. In particular, one can translate each line above from English to Haskell and get a working example. Let see.

  • A contestant in a game show faces three doors.
    data Door = Door1 | Door1 ">Door2 | Door3 deriving (Show, Eq, Enum)
    doors = [Door1 .. Door3]
    
  • There is a car behind one of them and goats behind the other two.
    data Object = Car | Goat deriving (Show, Eq)
    type Universe = [Object]
    
    nature :: MC Universe
    nature = shuffle 3 [Car, Goat, Goat]
    

Here shuffle is a function from monte-carlo package, which as the name suggests shuffles the first three objects from [Car, Goat, Goat]. This represents our state of the world.

  • He picks a door
    contestant :: MC Door
    contestant = sample 3 doors

The function sample is also from monte-carlo package. It picks one out of the list of options.

  • and then the game show host, who knows what is behind the doors, opens one of the other doors to reveal a goat.

Lets break down each English phrase. First, we need to figure out what options the game show host has. He can open any door which has a goat behind it. So, first we figure out which of the remaining doors have a goat behind them.

 options :: Universe -> Door -> [Door]
 options universe door = filter (\d -> Goat == open universe d) remainder
   where remainder = delete door doors

 open :: Universe -> Door -> Object
 open [a,_,_] Door1 = a
 open [_,a,_] Door2 = a
 open [_,_,a] Door3 = a
 open _        _    = error "Wrong number of elements in the universe"

Now the host can pick any one of the them.

 host :: [Door]-> MC Door
 host o = sample (length o) o where
  • The contestant then has the option to stick to his original choice, or switch to the unopened door.

I will write a function for each strategy of the contestant, even though one function is redundant.

 strategyStick :: Door -> Door -> Door
 strategyStick d1 ">d2 ">d1 d2 = d1
strategyFlip :: Door -> Door -> Door
 strategyFlip d1 d2 = head (doors \\ [d1,d2])
  • What is better?

Ah, the finale. Assuming that the contestant wants a car over a goat, we need to associate a reward with the outcome. If the outcome is a Car, the contest gets a reward of 1 Car, otherwise he gets 0 Cars.

 reward :: Universe -> Door -> Double
 reward u d = if Car == open u d then 1.0 else 0.0

Now, which strategy gives the contestant a better reward on average. We can, of course evaluate the average reward of each strategy analytically. We can also use the probability monad by Martin Erwig, which does exact probability calculations in Haskell. It even has Monty Hall problem as one of the examples. In this post, however, I am interested in a Monte Carlo approach. Monte Carlo approach works on the law of large numbers: one can approximate the expected value of a random variable by taking independent samples from its distribution and computing the sample average. To do so in Haskell, I first need to define the random variable, which in this case is the reward that we get. We need to specify the strategy in order to determine the outcome.

 randomVariable :: (Door -> Door -> Door) -> MC Double
randomVariable strategy = do
   u <- nature
   c <- contestant
   h <- host (options u c)
   let r = reward u (strategy c h)
   return r

This function says that the outcome is produced as follows. First nature chooses the ordering of the objects behind the door, then the contestant chooses a door, then the host opens a door, and then, depending on the contestant’s strategy, he gets a reward.

This is another good thing about Haskell. I can switch to an imperative style whenever it is convenient. In this case I can sequence a bunch of random events, and I have not even talked about generating random number!. Haskell provides a complete separation of generating random variable from the Monte Carlo simulation. Thanks to the new repeatMC function in the monte-carlo package, running a simulation is straight-forward.

 main = let
   n     = 100000
   seed  = 42
   stats = repeatMC n (randomVariable strategyStick) `evalMC` mt19937 seed
   in do
putStrLn $ printf "Mean                    : %f" (sampleMean stats)
      putStrLn $ printf "99%% Confidence interval : (%f, %f)"
                `uncurry` (sampleCI 0.99 stats)

The results of the simulation are stored in the Summary data type, which I talked about in the last post.

I only output the result of the stick strategy. The result of the switch strategy is going to be (1 - result of stick), since if the contestant did not get the car by sticking to his choice, he would have gotten it by switching. And the result (surprise, surprise, if you did not know about Monty Hall problem)

Mean                    : 0.3355999999999993
99% Confidence interval : (0.33175368331333815, 0.3394463166866604)

So, sticking to one’s choice only wins 1/3rd of the time. Which means that switching wins 2/3rd of the time. So, if you happen to be in a Monty Hall game show—switch.

Summary data type in monte carlo monad

The new version of Haskell package monte-carlo by Patrick Perry adds a repeatMC function which is useful for some Monte Carlo experiments. In the course of implementing this, Patrick has also implemented a Summary data type for keeping track of basic statistical properties of Monte Carlo experiments. This data type is also useful on its own, though there is no interface to use it independently. In this post I show how it can be used outside the MC monad.

import Control.Monad.MC
import Data.List (foldl')
import Text.Printf (printf)

Suppose we have some data in a list of Doubles that we want to analyze.

x = [1..10] :: [Double]

First we need to create a summary from the data. There is no pre-defined function for it, but we can easily define one on our own.

stats :: [Double] -> Summary
stats = foldl' update summary

Now, we can use this functions to get the basic statistics of the data.

s = stats x

There are a few functions to extract the statistics from the summary. See Control.Monad documentation for the complete list and description. Most useful ones for me are those that give basic properties of the data

  • sampleSize :: Summary -> Int
  • sampleMin :: Summary -> Double
  • sampleMax :: Summary -> Double

and those that give basic statistics of the data.

  • sampleMean :: Summary -> Double
  • sampleVar :: Summary -> Double
  • sampleSD :: Summary -> Double

We can obtain these functions as follows:

main = do
putStrLn $ printf "Data   : %s" (show x)
  putStrLn $ printf "Length : %d" (sampleSize s)
  putStrLn $ printf "Max    : %f" (sampleMax s)
  putStrLn $ printf "Min    : %f" (sampleMin s)
  putStrLn $ printf "Mean   : %f" (sampleMean s)
  putStrLn $ printf "Var    : %f" (sampleVar s)
  putStrLn $ printf "SD     : %f" (sampleSD s)

which gives

Data   : [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
Length : 10
Max    : 10.0
Min    : 1.0
Mean   : 5.5
Var    : 9.166666666666666
SD     : 3.0276503540974917

Sometime next week I will show how to use this data type in Monte Carlo simulations.