I want to begin our coding adventure with something simple yet useful in many areas. I think that has to be using Monte Carlo simulation to solve a problem, especially that for which we do not know the exact solution.

In essence, there are 2 stages for implementing a Monte Carlo simulation:

1. Input stage: Define a set of possible inputs for the problem. Randomly generate inputs from this set. We can do this one at a time, or altogether.
2. Output stage: For each generated input, find the output by following the problem description. I will refer to this “give me an input and I’ll give you an output” as a simulation. Combine the outputs to arrive at some meaningful result.

The benefit of a computer program shines in both stages. Our computer can create many inputs and evaluate their outputs much faster than we can.

Over the next few posts, we will examine how we can apply Monte Carlo simulation to solve problems in probability and integration. Whenever possible, we will consider the exact solution and see how good our approximate solution is.

1. Problem description

In a game of craps, we roll two dice and find their sum. There are three possible scenarios leading to a win or loss.

1. If the sum is 7 or 11 (called natural), we win.
2. If the sum is 2, 3, or 12 (called craps), we lose.
3. Otherwise, we roll the dice until we get the initial sum (win) or a sum of 7 (lose).

What is the probability that we win a game of craps?

Before we start coding, let us first agree upon the input and the output of a simulation. We are interested in the probability of winning craps. Finding the exact value (i.e. with a proof) seems difficult. However, it seems reasonable that if we play craps times, then the ratio

will be similar to the exact probability of winning craps.

is a fixed number, e.g. 100,000, but the number of wins can change. Therefore, the output of a simulation should be whether we won or lost a game of craps. By doing so, we can tally how many times we won.

What is the input then? At first, it seems like the input is the sum of the dice. I’d argue that we should be more rigorous, however, because there may be times in which we need to throw additional rolls (scenario 3). Let us define the input to be the sequence of sums. Then, sequences such as , , , and will result in a win, and sequences like , , , and will result in a loss. Perfect.

2. A single game

I really like this problem because we can directly translate the rules of the game to a code, just using basic programming concepts. The code below shows a set of Matlab commands to run a game of craps:

% Roll the dice and find their sum
die1 = randi(6);
die2 = randi(6);
sum  = die1 + die2;

% Check if we have a natural
if (sum == 7 || sum == 11)
    display('We won!');

% Check if we have craps
elseif (sum == 2 || sum == 3 || sum == 12)
    display('We lost...');

% Otherwise, we roll the dice until we get the initial sum
% or a sum of 7
else
    while (true)
        % Roll the dice and find their sum
        die1    = randi(6);
        die2    = randi(6);
        sum_new = die1 + die2;

        if (sum_new == sum)
            display('We won!');

            break;

        elseif (sum_new == 7)
            display('We lost...');

            break;

        end
    end
end

Notice how the if-else statements on lines 7, 11, 16, 23, and 28 accurately reflect the win and loss conditions. Line 17 shows a while loop. While loops are great when we do not know how many times we need to repeat something. In particular, the statement “while (true)” allows us to hide the details of (i.e. not worry about) when the loop should stop running. We can do that later with a break statement, after we have all the big pieces in place.

One thing that I glossed over is the input of a simulation. Before, we had defined the input to be the sequence of sums. However, nowhere in our code above do we find that sequence. At least, not in its entirety. Going back to the rules of craps, we see that we do not need to remember all numbers in the sequence, but just need to keep track of the initial sum (i.e. the first number in the sequence) and the new sum that will end the simulation (i.e. the last number in the sequence).

3. Monte Carlo simulation

We can play a single game of craps. To estimate the probability of winning craps using Monte Carlo simulation, we need to play the game multiple times, say times. Recall that a loop is great for repeating something. Since we know how many times, let’s use a for loop.

Here is a pseudocode in Matlab:

% Reset the number of wins
numWins = 0;

% Set the number of simulations
N = 10^5;

% Run the simulations
for i = 1 : N
    % Run a single game of craps and increase the number
    % of wins by 1 if we win the game
    [...]
end

% Display the probability of winning craps
display(numWins / N);

The probability that I got is 0.49228, slightly less than 0.5. It seems that the game is unfavorable to the player.

4. Mathematical proof

Let’s find the probability of winning craps. There are two possible ways to win the game: (1) from the initial roll, when the sum is 7 or 11, or (2) from a subsequent roll, when the new sum equals the initial sum.

Hence,

For convenience, I denoted the two probabilities on the right-hand side as (for the first case) and (for the second case).

Calculating the probability that the first case occurs is easy:

But what about the second case, when we can’t know how many rolls are needed to get the initial sum again? It’s not hard, it turns out.

Suppose the initial sum is 4. Let denote the probability that, in a subsequent roll, we get a sum of 4 before we get 7. For brevity, let’s omit the conditional clause “initial sum is 4.” Can you see why the following equation is true?

There is a chance that we get a 4 on the first try (technically, second roll). If we did not get a 4, then as long as we did not get 7 either (this occurs with the probability of ), we get another try at rolling the dice.

Solve for the variable to get,

Hence, the probability that we win from an initial sum of 4 is,

Similarly, we can compute the probabilities that we win from an initial sum of 5, 6, 8, 9, or 10. Add all these to get :

We conclude with the probability of winning craps:

Compare my approximate value of 0.49228 to this exact value. There is only 0.132% relative error with my guess! By running more simulations, we can make the relative error even smaller.

Notes

I was inspired to write a Monte Carlo simulation code for craps when I TAed for an introductory programming class. I thought it made a great problem for the students after they learned to use a while loop. (It’s not easy to come up with a problem that can be solved with a while loop but can’t/shouldn’t be solved with a for loop!)

John Haigh wrote a book called “Taking Chances,” which focuses on probability in games. He derives the probability that I called using a different argument:

“If events A and B are disjoint, so they cannot occur simultaneously, and their respective chances are a and b, then the chance that A comes first is a / (a + b).

To see why that is so, note that if neither A nor B occurs, we can simply ignore that throw. The only throws that are relevant are those that result in either A or B. Because their respective probabilities are a and b, in a long series of N throws A and B occur about Na and Nb times. Concentrating only on those throws that give either A or B, the occurrences of A and B are mixed up at random; the chance that A comes first is just the proportion of As in this list—and that is Na / (Na + Nb) = a / (a + b).”

You can find the code in its entirety and in other languages here:
Download from GitHub