The Mathematics of Modeling
Part 2: Monte Carlo Methods

Our goal in computational science is to build a model that realistically reflects some particular behavior of an actual, physical system.  Occasionally the factors that affect the behavior are so numerous and varied that it is easier to consider the possibilities as the results of random chance, favored by particular probabilities.  Problems that involve a random variable are called stochastic problems, and modeling methods that take advantage of random variables are called Monte Carlo methods, in reference to the famous casino and resort town in Monaco.

LCG Random Number Generation
	The key element to a Monte Carlo method is the ability to generate a random number.  There is actually quite a bit of philosophical discussion about what makes a number truly 'random' - computers are naturally deterministic, i.e. what comes out is pre-determined by what goes in, making it quite a challenge for programmers to create 'pseudo-random' algorithms to generate sequences of numbers that sufficiently fit mathematical definitions of randomness. Let's take a look at one common method, called a Linear Congruential Generator (LCG).

An LCG algorithm requires 4 starting numbers.  The first number is known as the seed, because it acts as the starting point for a list of numbers, and only gets used to generate the first number. Computers will often choose this number from some other sort of 'random' source (typically a time from the computer's clock) to add even more variability.

seed  =  17;

The other three numbers - modulus, scale and shift - are used for each step.

scale  =  5;
shift  =  2;
modulus  =  13;

We multiply the seed by the scale...

num  =  seed*scale

Add the shift...

num  =  num  +  shift

And take the remainder after dividing by the modulus. This gives us a number that is greater than or equal to zero, but less that the modulus. Mathematica uses the function Mod for modulo (remainder) division.

num  =  Mod[num, modulus]

The new number then becomes the seed for the next random number...

seed  =  num;
num  =  Mod[scale*seed + shift, modulus]

We can make a function that takes our four numbers as arguments (as well as a list length) and gives us a list of pseudorandom numbers between (and including) zero and (but not including) the modulus.

LCG[length_, seed_, scale_, shift_, 
        modulus_]  := 
(num  =  
          seed; 
Table[
num  =   Mod[
              scale*num + shift, modulus]
   , {i, 
            length}])

LCG[20, 17, 5, 2, 13]

Not very random, is it! Let's try a different scale...

LCG[20, 17, 7, 2, 13]

Yikes! that's even worse! But if we keep the 'bad' scale and change the shift...

LCG[20, 17, 7, 3, 13]

Which initially looks better, until you realize that the series again starts to repeat, after selecting every number from 0 to 12 (except for 6 !?!).  A major problem with this kind of generator is that once a number is repeated in a series, the entire series begins to repeat!  We can make the distance between 'repeats' longer with a bigger modulus...

LCG[20, 17, 7, 3, 17]

But the patterns still appear, unless we decide not to use just integers...

LCG[20, 17.2, 7.3, 3.01, 17.4]

Now these numbers seem much more random, at least at first glance. So the question is, are the patterns 'gone' or are they just harder to see? One way to get a picture of what is happening is to plot each random number along with one of it's neighbors on an x-y graph! We can do this by Partitioning the list into groups of 2 and then using ListPlot.

RandList  =  
        Partition[ LCG[20, 17.2, 7.3, 3.01, 17.4], 
          2];
ListPlot[RandList, Prolog -> AbsolutePointSize[4]]

The points appear to be fairly spread out ('uniformly distributed'),  without an apparent pattern - what if we take more points?

RandList  =  
        Partition[ LCG[200, 17.2, 7.3, 3.01, 17.4], 
          2];
ListPlot[RandList, Prolog -> AbsolutePointSize[4]]

Apparently there is still an underlying pattern (see why it is called a Linear Congruential generator?), however, another choice of numbers looks much better for 200 numbers...

RandList  =  
        Partition[ LCG[200, 17.2, 45.3, 3.01, 17.3], 
          2];
ListPlot[RandList, Prolog -> AbsolutePointSize[4]]

But patterns still start to appear if we take enough points!

RandList  =  
        Partition[ LCG[2000, 17.2, 45.3, 3.01, 17.3], 
          2];
ListPlot[RandList, Prolog -> AbsolutePointSize[4]]

Perhaps you can see now why there is alot of study and discussion about how to generate random numbers - even the most complex algorithms may well have underlying patterns that belie randomness.

Mathematica has it's own 'uniformly distributed pseudorandom' number generator, contained in the function Random.  To simulate rolling a pair of dice, we ask Mathematica to generate random integers from the range of 1 to 6.  We use the := notation to tell Mathematica to 'evaluate later', that way, each time we ask for a DieRoll, Mathematica will evaluate the Random function and select a new random integer.

DieRoll  :=  Random[Integer, {1, 6}]

{DieRoll, DieRoll}

In the casino spirit, we can also simulate drawing one card from a deck of cards by randomly choosing one of the four suits and one of the thirteen ranks.  The integer we choose represents the position of the rank or suit in the appropriate list.

\(\(Suits\  = \ {"Clubs", "Diamonds", "Hearts", "Spades"};\)\)
\(\(Ranks\  = \ {"Ace", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine", "Ten", "Jack", "Queen", "King"};\)\)
Suitpick  :=  Random[Integer, {1, 4}]
Rankpick  :=  Random[Integer, {1, 13}]

Card  =  
      Ranks[[Rankpick]]  <> " of " <> 
        Suits[[Suitpick]]

Estimation of Pi
	The 'classic' problem for introducing Monte Carlo methods is a model for estimating the value of π. The model is meant to illuminate one simple concept - that the relationship between the area of a circle and the area of a square is indeed subject to this 'mysterious' factor of approximately 3.14159...
	According to geometry, the area of a circle is subject to the formula, (Subscript[A, circle])  =  ((π (d^2))/4), where d is the diameter of the circle.  If the circle is inscribed within a square, the sides of the square will be length d, and the area of the square (Subscript[A, square]) =  (d^2). The ratio of the two areas is theoretically ((Subscript[A, circle])/(Subscript[A, square])) =  (π/4).

Show[
      Graphics[
{Rectangle[{-1, -1}, {1, 
              1}], 
GrayLevel[
            0.7], 
Disk[{0, 0}, 
            1], 
GrayLevel[
            0.0], 
Line[{{-1, 0}, {1, 
                0}}]
}], 
      AspectRatio  ->  Automatic]

To use a Monte Carlo method, we use the picture above as a sort of 'dartboard' and randomly select points in the entire square, some of which will also lie inside the circle. We choose these points as real numbers from -radius to +radius, thus setting the diameter of our circle.

Diameter  =  2;  Radius  =  Diameter/2;
Ranx  :=  
        Random[Real, {-Radius, Radius}];
Rany  :=  
        Random[Real, {-Radius, Radius}];
{Ranx, Rany}

We now make a list of random points using the Table function, and look at the first few elements using Take.

AllPoints  =  500;
PointList  =  
        Table[{Ranx, Rany}, {i, AllPoints}];
Take[PointList, 5]

What we care about for each point, though, is how far away it is from the origin.  Points that are closer to the origin than the radius of the circle are in the circle, points that are not are outside the circle.  We use a For loop to go through each point (from 1 up to our variable AllPoints) and count the number in the circle and out of the circle, using the distance formula as the criterion.

InCirc  =  0;
OutCirc  =  0;
For[i = 1, 
        i <= AllPoints, i++, 
{
xval  =  PointList[[i]][[1]]; 
yval  =  PointList[[i]][[2]]; 
Dist  =  (Sqrt[(xval^2) +  (yval^2)]); 
If[
            Dist < Radius, InCirc++, OutCirc++];
}];
InCirc
OutCirc

If we we plot each point, we can see what our results look like, although it is probably a good thing that we had the program count which points were in and which were out...

Show[ ListPlot[PointList], 
        Graphics[
          Circle[{0, 0}, Radius]], 
AspectRatio -> 
          Automatic];

Probability dictates that the proportion of points that lie in the circle will correspond to the ratio of areas, i.e. ((Subscript[p, in])/(Subscript[p, total])) = ((Subscript[A, circle])/(Subscript[A, square])) =  (π/4). Given that we know (Subscript[p, in])and (Subscript[p, total]), we can make an estimate for π, where \(TraditionalForm\`π\_est = \(\(4\) \( p\_in\/p\_total\)\(\ \)\)\).

PiGuess  =  N[(4*InCirc)/AllPoints]

Not a bad estimate!  Probability tends to work even better 'in the long run', so we would expect better results if we select more random points.  This is much easier if we compile everything we've done into one guessing function that takes the number of points as an argument, and returns an estimate:

\(PiGuessFn[NumPoints_]\  := 
\((Radius\  = \ 1.0; 

        InCirc\  = \ 0; 
OutCirc\  = \ 0; 

        PointList\  = \ 
          Table[
\((Ranx\  = \ 
                Random[Real, {\(-Radius\), 
                    Radius}]; 
Rany\  = \ 
                Random[Real, {\(-Radius\), 
                    Radius}]; 
Dist\  = \ \@\(Ranx\^2 + \ Rany\^2\); 
If[
                Dist < 1, \(InCirc++\), \(OutCirc++\)]; 
{Ranx, Rany})\), {NumPoints}]; 

        N[\((4*InCirc)\)/NumPoints])\)\)

The general trade-off in computational modeling is that getting better results takes more computer time. Calculating with one hundred thousand points might take a bit longer, but should give us an even better estimate!

PiGuessFn[100000]

N[Pi]

The trade-off between accuracy and expense is interesting to many computational scientists.  Because we have a good 'correct' answer to the Pi estimation, we can determine how much error each of our simulations has.  We can then demonstrate how a sequence of increasingly 'expensive' (in terms of computer time) simulations will follow a general trend of smaller errors.

PiGuessList  =  
      Table[{2^i, Abs[Pi - PiGuessFn[2^i]]}, {i, 2, 18}]

ListPlot[PiGuessList[[All, 2]], 
      PlotJoined -> 
        True, 
AxesLabel ->  {"Expense", "Error\
"}]

This notebook created by Dr. Garrett R. Love for the Shodor Education Foundation. Copyright 2001