Using Simulations to Answer More Questions

As the modules above execute, a lot more information can be gathered. For example, in the second module we kept track of the number of times the gambler played the game, but did nothing with that information. Suppose that we are asked to determine the average number of times a gambler will play if they win or lose. These values would be important to the casino in terms of how fast they will make money. In this section we will put together a more complicated Module that will perform the simulation as in the previous section and return lists of numbers that can be used to show the distribution of the number of times the gamblers played in the simulator.

One way to keep track of the gamblers is in terms of the distribution of the number of times winning or losing gamblers play. That is, different winning gamblers will play a different number of time. We can keep track of the number of gamblers that win in exactly 10, 11, 12, etc. plays of the game. In the same way we can keep track of the distribution of losing gamblers in terms of the number of times they play before being ruined. The module that must tally up the number of gamblers in each category. An example of a module that will compute and return these lists is the following.

  SGR[ a_, c_, p_, ngamblers_, maxnum_ ] :=
     Module[ {output,nplays,flist,wlist},
           wlist = {};
           flist = {};
           Do[AppendTo[wlist,0], {i,1,maxnum}];
           Do[AppendTo[flist,0], {i,1,maxnum}];
           Do[Clear[output,nplays];
             output = GamblersRuin[a,c,p];
             nplays = output[[1]];
             If[nplays <= maxnum,
                If[output[[2]],
                   wlist[[nplays]] = wlist[[nplays]] + 1,
                   wlist[[nplays]] = wlist[[nplays]] + 1
                   ]
               ],
             {i,1,ngamblers}
             ];
             Return[{tlist,flist}]
             ]

The module is built by combining the ideas in the first modules developed and then including a little more code to keep track of the results. Note that we again are using the original module, GamblersRuin, to do the basic simulation of a single gambler. The features in the new module are the following.

1. Another input parameter has been added. The variable is used to keep track of the longest allowable sequence. Since a gambler could be in the casino for thousands of plays before winning or losing, this variable ignores any sequences that are too long. This is done for practical reasons. We cannot store arbitrarily long sequences due to finite resources (memory) on a computer.
2. There are two list variables, wlist and flist, that are used to keep track of the numbers. Note that the two lines

              wlist = {};
              flist = {};
   

   just initialize the two variables to empty lists.
3. The next two lines,

              Do[AppendTo[wlist,0], {i,1,maxnum}];
              Do[AppendTo[flist,0], {i,1,maxnum}];
   

   create and initialize the lists to zeros. Thus before the simulation begins, there are two lists of length maxnum with zeros in all the entries. If the gambler wins in exactly seven plays, then wlist[[7]] will be incremented by +1.
4. The Do command then loops over all the gamblers and does the simulation. Notice that in this case, the number of times the gambler plays is used. The first value returned by the GamblersRuin module is the number of plays. If the number of times played is larger then the maximum sequence length, maxnum, nothing is done to either list. If the value is in range the appropriate location in the list is incremented by one to indicate that another gambler has won or been ruined in a given number of plays. The location in the list corresponds to the number of plays.

The code is a little more complicated, but the pieces should make sense. The main reason for limiting the number of times an individual gambler can play is for computational reasons. If we allowed the gambler to play for longer than some finite number, the simulator would take a long time to run and the size of the lists may become very large. In fact, it can be shown that the probability of a very long sequence is small. That is, most of the sequences that gamblers would encounter are relatively short. Note that this is a place where we are again making an approximation that must be justified. We can use the simulator to help us with this; see Problem 9.

To try the new module out we can use the command

  simlist = SGR[5, 10, 0.44, 50, 40]

The input indicates that the gamblers are staring with $5, are trying to end with $10, and are playing a game with a probability of 0.44 of winning. The command above will do the simulation for 50 gamblers and will throw out any cases where the gambler plays more than 40 times. The output from the module is a pair of lists. The first list contains the data from the cases where the gambler wins and the second list contains the information for the cases when the gambler loses. You should note that if you execute the same command above your results will be different since Mathematica produces a different random sequence each time the routine is executed.

Figure 1.2: Barcharts from the simulation of the Gambler's Ruin problem showing the distribution of the number of times the gambler wins and loses in terms of the number of times the game is played. The third barchart shows both distributions together. This simulation was done with 50 computer simulated gamblers.

We can visualize the results above by creating a histogram of the data. To make histograms we need to extract each of the lists and then use the graphics on the separate lists. To do this we can use the Part command in Mathematica. We could use

  winlist = Part[simlist,1]

to grab the list describing the distribution of the number of plays it takes gamblers to win. The output will look something like

  Out[6]= {0, 0, 0, 0, 4, 0, 6, 0, 7, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 3, 
  >    0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 2, 0, 0, 0, 1, 0}

The list defining the distribution of the losing simulated gamblers can be obtained using the command

  loselist = Part[simlist,2]

Now that we have both lists in the separate variables, winlist and loselist, we can work with the two sets of data.

To create a histogram, we can use the BarChart command as described in Appendix (). Note that the BarChart command is part of a package that must be loaded. You can do this via the function/package browser under the Help menu on the notebook (load the Graphics package). To get an appropriate scaling for the histograms, we need to get the maximum values of each of the lists. These can be grabbed and stored into variables using a command like

   maxval = Max[winlist]

These values will be used to set the PlotRange for the graphics we will do. The commands to create the histograms are

 maxval = Max[winlist]
 BarChart[winlist, PlotRange -> {0,maxval}]

and

 maxval = Max[loselist]
 BarChart[loselist, PlotRange -> {0,maxval}]

Typical examples of the barcharts are shown in Figure (1.2).

Figure 1.3: Barcharts from the simulation of the Gambler's Ruin problem showing the distribution of the number of times the gambler wins and loses in terms of the number of times the game is played. The third barchart shows both distributions together. This simulation was done with 500 computer simulated gamblers.

Figure 1.4: Barcharts from the simulation of the Gambler's Ruin problem showing the distribution of the number of times the gambler wins and loses in terms of the number of times the game is played. The third barchart shows both distributions together. This simulation was done with 5000 computer simulated gamblers.

If we performed a simulation that simulates 500 gamblers using

  sim500 = SimulateGR[5,10,.44, 500, 40]

we should obtain better estimates of the probability and other information about the process. The barcharts for this case are shown in Figure (1.3). From these barcharts we are getting a better reflection of the distribution of events. It is also easy to see that the probability of winning is not very high in this case. That is, there appear to be a lot more losers than winners. To do an even better job we might consider running the simulator for 5000 gamblers using

  sim5000 = SimulateGR[5,10,.44, 5000 ,40]

These barcharts are shown in Figure (1.4). Again, we see an improvement in the results. From these histograms we should start getting an idea of how the number of times played by a gambler is distributed.