Monte Carlo Simulation

In order to make statistical predictions about the long-term results of a random process, it is often useful to do a simulation based on one's understanding of the underlying probabilities.  This procedure is referred to as the Monte Carlo method.

As an example, consider a casino game in which a player bets against the house and the house wins 51% of the time.  The question is, how many games have to be played before the house is reasonably sure of coming out ahead.  This scenario is common enough that mathematicians long ago figured out very precisely what the statistics are, but here we want to illustrate how to get a good idea of what can happen in practice without having to absorb a lot of mathematics.

First we define an expression that computes the net revenue to the house for a single game, based on a random number chosen from 1 to 100.  If the random number is 1 to 51, the house wins one betting unit, whereas if the number exceeds 51, the house loses one unit.  (In a high-stakes game, each bet may be worth $1000 or more.  Thus it is important for the casino to know how bad a losing streak it may have to weather in order to turn a profit--so that it doesn't go bankrupt first!)

revenue  :=  
      If[Random[Integer,  {1,  100}]  <=  51,  1,  -1]

We can play the game by evaluating the expression revenue.  The following command plays the game 10 times and lists the results; i is a dummy variable that runs from 1 to 10.

Table[revenue,  {i,  10}]

The net profit to the casino after 10 games is found by adding the 10 numbers on the output line above.  Next we create a function that gives the net profit for the house after n games.

profit[n_]  :=  Sum[revenue,  {i,  n}]

On average, every 100 games the house should win 51 times and the player(s) should win 49 times, so the net profit to the house should be 2 betting units.  Let's see what happens in a few trial runs.

Table[profit[100],  {i,  10}]

We see that the net profit can fluctuate significantly from one set of 100 games to the next, and there is a sizable probability that the house has lost money after 100 games.  To get an idea of how the net profit is likely to be distributed in general, we can repeat the experiment a large number of times and make a histogram of the results.

<< Statistics`<< Graphics`

Having loaded the Statistics and Graphics packages, we write a function that plays the game n times, records the net profit, and repeats the experiment k times.  It then tabulates the results, setting up intervals, or "bins", of length xstep centered at x0, x0 + xstep, x0 + 2xstep, ..., x1, counting how many times the net profit after n games falls into each interval, and graphing the result as a histogram, or bar chart.  The idea is to figure out, in a systematic way, what possible outcomes the casino can reasonably expect after n games.

hist[n_,  k_,  x0_,  x1_,  xstep_]  :=  
      GeneralizedBarChart[
        Transpose[{Table[x,  {x,  x0,  x1,  xstep}],  
            BinCounts[Table[profit[n],  {i,  k}], 
               {x0  -  xstep/2,  x1  +  xstep/2,  xstep}],  
            Table[xstep,  {x,  x0,  x1,  xstep}]}],  
        PlotRange  ->  All,  AxesOrigin  ->  {x0  -  xstep/2,  0}]

Now we run the function hist using n = 100 and k = 100.  Theoretically the house could win or lose up to 100 chips, but in practice we find that nearly all of the outcomes are within 30 (on either side) of the expected value of 2.  We choose the upper and lower limits of our histogram accordingly.  The vertical axis below represents the number of times a given value of the net profit is obtained during the 100 trials.

hist[100,  100,  -28,  32,  2];

The histogram confirms our impression that there is a wide variation in the outcomes after 100 games.  It looks like the casino is about as likely to have lost money as to have profited.  However, the distribution shown above is irregular enough to indicate that we really should run more trials to see a better approximation to the actual distribution.  Let's try 1000 trials.

hist[100,  1000,  -28,  32,  2];

According to the Central Limit Theorem, when both n and k are large, the histogram should be shaped like a "bell curve", and we begin to see this shape emerging above.  Though it will take a while to run, let's move on to 10,000 trials.

hist[100,  10000,  -28,  32,  2];

Here we see very clearly the shape of a bell curve.  Though we haven't gained that much in terms of knowing how likely the house is to be behind after 100 games, and how large its net loss is likely to be in that case, we do gain confidence that our results after 1000 trials are a good depiction of the distribution of possible outcomes.

Now we consider the net profit after 1000 games.  We expect on average the house to win 510 games and the player(s) to win 490, for a net profit of 20 chips.  Again we start with just 100 trials.  Since both the range of possible outcomes and the expected outcome have increased by a factor of 10, we may be inclined to increase the horizontal range of the histogram by a factor of 10 as well.

hist[1000,  100,  -280,  320,  20];

We find that, relatively speaking, the outcomes are clustered much closer together.  This reflects the theoretical principle (also a consequence of the Central Limit Theorem) that the average "spread" of outcomes after a large number of trials should be proportional to the square root of n, the number of games played in each trial.  This is important for the casino, since if the spread were proportional to n, then the casino could never be too sure of making a profit.  When we increase n by a factor of 10, the spread should only increase by a factor of (Sqrt[10]), or a little more than 3.  Let's adjust the range of our histograms accordingly.

hist[1000,  100,  -80,  120,  8];

We see that after 1000 games, the house is definitely more likely to be ahead than behind.  However, the chances of being behind are still sizable.  Let's repeat with 1000 trials to be more sure of our results.

hist[1000,  1000,  -80,  120,  8];

We see the bell curve shape emerging again.  Though it is unlikely, the chances are not insignificant that the house is behind by more than 50 chips after 1000 games.  If each chip is worth $1000, then we might advise the casino to have at least $100,000 cash on hand to be prepared for this possibility.  Maybe even that is not enough--to see we would have to experiment further.

Finally, let's see what happens after 10,000 games.  We expect on average the house to be ahead by 200 chips at this point, and based on our earlier discussion the range of values we use to make the histogram need only go up by a factor of 3 or so from the previous case.  Even 100 trials will take a while to run now, but we have to start somewhere.

hist[10000,  100,  -100,  500,  20];

It seems that turning a profit after 10,000 games is highly likely, though with only 100 trials we do not get such a good idea of the worst-case scenario.  Though it will take a good bit of time, we should certainly do 1000 trials or more if we are considering putting our money into such a venture.  If you run the following command, you may want to take a long break from the computer and come back later to see what happened.

hist[10000,  1000,  -100,  500,  20];

Though the chances of a loss after 10,000 games is quite small, the possibility cannot be ignored, and we might judge that the house should not rule out being behind at some point by 100 or more chips.  However, the overall upward trend seems clear, and we may expect that after 100,000 games the casino is overwhelmingly likely to have made a profit.  Based on our previous observations of the growth of the spread of outcomes, we expect that most of the time the net profit will be within 1000 of the expected value of 2000.  We show the results of 10 trials of 100,000 games below.

Table[profit[100000],  {i,  10}]