Studying craps using Haskell

from Data.Random by Alex Mason

Last week, a friend and I were working on a piece of homework for our engineering systems analysis course. We were supposed to be studying the probabilities of a player winning a game of craps, and then studying the effect of using the martingail betting system (double your bet if you lose, start again with your starting bet when you win).

We were supposed to be using MATLAB for this, but neither of us for the life of us could figure out how we were supposed to be able to encode that system in MATLAB. We spend an hour or so trying to figure it out, and then I realised that, while we could use an array/matrix, the game was perfectly modelled using a binary tree.

We went looking through the docs to see what support for trees MATLAB had; it seemed to have some, but neither of us could figure out how to use them (I am really starting to dislike MATLAB's docs, there's usually one sentence on a page that's useful, and it's very hard to get all the info you need... but that's another matter). I've done quite a bit of work using binary trees in haskell, so I decided I'd try and get things working in Haskell.

I started out with this datatype:

data Tree = Node Integer Tree Tree | Leaf Integer
            deriving (Show)

where each Node would hold the current profit, the left would be the tree for losing, and the right for winning. The Leaves hold the end profit after n games. To create the tree, I created a buildTree function:

buildTree 0 p = Leaf p
buildTree n p = Node p (buildTree (n-1) (p*2)) (buildTree (n-1) p)

Hooray for bugs, that's totally not what the rules said we're supposed to do (it wasn't until a few hours later we found out we had interpreted the rules incorrectly, so this definition was correct according to out version of the rules). Anyway, you can see that creating a tree of profits is pretty damn easy. We also had to keep track of the amount of money bet to get to that result, and calculate the probability of reaching that result. Since we're working with a purely Bernoulli distribution (probability of winning to the power of the number of wins times the probability of losing (1-prWin) to the power of the the total number of games minus the number of wins), I made a function to calculate that:

binomial prw num wins = prw^(wins) * (1-prw)^(num-wins)

To hold all this info we needed, I created a new datatype

type Profit = Int64       -- Profits are large integers
type Stake = Int64        -- So are stakes
type Probability = Double -- Probabilities are floating point numbers
type Bets = Int64 
data P = P Profit Stake Probability
         deriving (Eq,Ord)

And once I'd done that, I needed to rewrite my buildTree function to actually produce the correct values. In this function, I used a helper function, which would keep track of the the net profits, the total stake bet, and compute the probability of this happening based on the number of times we won getting to that point in the tree (prWin was calculated in our lab session using MATLAB):

prWin, prLoss :: Double
prWin = 0.492929292929293
prLoss = 1 - prWin

-- buildTree forms the tree of all wins and losses to a given
-- depth n along with the profit after n games, and the total
-- stake paid during the game.
buildTree :: Bets -> Stake -> Tree
buildTree num stake =
        buildTree' num 0 0 stake stake
        -- buildTree' is the helper function that actuall creates the tree
        where -- when n == 0, we have played n games, and we put the result
              -- leaves in the tree
              buildTree' 0  prof wins _ statot
                = Leaf (P prof (statot-stake) (binomial prWin num wins))
              -- In n != 0, then we create a Node with the tree after
              -- losing on the left, and the tree after winning on the
              -- right.
              buildTree' n' prof wins sta statot
                = Node (buildTree' (n'-1) (prof - sta) wins     (sta * 2) (statot+sta))
                       (buildTree' (n'-1) (prof + sta) (wins+1)  stake    (statot+sta))