Estimate Pi Using Monte Carlo Sampling

A classic use of the Monte Carlo random sampling technique is to estimate Pi. It's very simple.

1. We construct a square by randomly selecting point in the interval (-1, 1).

2. We select a subset of points that satisfy the predicate that places them on or within the unit circle:

x^2 + y^2 < =1.

3. The area of the square extending from -1 to 1 is 2^2 = 4. The area of the unit circle is 1 (Pi*r^2 = Pi*1^2). Therefore we multiply the ratio of the number of points lying within the unit circle divided by the number of points within the square by 4.

Here is an implementation. Below that is an implementation I found on the web that shows how newbies to Mathematica may transfer inefficient coding techniques from other languages to Mathematica.

Clear@circleRegionMonteCarlo;circleRegionMonteCarlo[numberSamplePoints_Integer]:=Module[{samplePoints=RandomReal[{-1,1},{numberSamplePoints,2}],
areaOfCirclePoints},

(* Select points within the unit circle using the parametric equation for a circle *)
areaOfCirclePoints=Select[samplePoints,First@#^2+Last@#^2<=1&];

Print@StringForm["Approximation to Pi using `` sample points is ``.",numberSamplePoints,4*Length@areaOfCirclePoints/numberSamplePoints//N];

If[numberSamplePoints<50001,ListPlot[areaOfCirclePoints,AspectRatio->Automatic]]
]

circleRegionMonteCarlo[50000]

Approximation to Pi using 50000 sample points is 3.14104`.

This represents an obsolete coding technique: 1) initializing the Pin (points inside the circle) and Pout (points outside the circle) arrays, 2) the For loop calling named random variable functions, 3) testing the distance predicate in the loop instead of all at once as is typical in functional programming, 4) having to take the Length of Pin and Pout in the loop (obviously you don't need to test the Length every time anyway), and 5) using Return instead of just leaving off the semi-colon after approx, are all obsolete techniques.

MonteCarloPi[n_]:=Module[{d,i},
n=n0;
Pin=Pout={};
For[i=1,i<=n,i++,
x=RandomReal[];
y=RandomReal[];
d=x^2+y^2;
If[d<=1, Pin=Append[Pin,{x,y}],Pout=Append[Pout,{x,y}]  ];
n=Length@Pin;
m=Length@Pout;
\[Rho]=m/n;
approx=4.0;]
Return@approx

]

Mathematica Plot - Basic Plotting

Plot

For starters, be aware that Mathematica's data plotting functionality is extensive and there are a variety of Plot functions. Using Information "?" with the wildcard characters "*" shows all functions that include "Plot". I won't show them here, but you can execute the command if you wish to see them.

?*Plot*

Basic Plotting

Let's say you have a function and you want to understand its behavior. We'll do a simple one.

In[238]:= y == x^n // TraditionalForm

Out[238]//TraditionalForm= y==x^n

Let's start by generating a few sample points of the function and plotting them.  You might think in the age of computers we can dispense with plotting sample points and just let Plot do the job. But when Plot fails to give you output and you don't know why, one measure is to go ahead and plot a few sample points by hand and compare them to what Plot is doing. So let's generate the sample points the easy way with Table. We Clear the function name, define the function, give x a value of 2, and use Table to sample the domain from -5 to 5 in increments of 0.5.

In[249]:= Clear@f; f[x_, t_] := x^t;

In[250]:= samplePoints = Table[f[2, t], {t, -5, 5, .5}]

Out[250]= {0.03125, 0.0441942, 0.0625, 0.0883883, 0.125, 0.176777, 0.25, 0.353553, 0.5, \
0.707107, 1., 1.41421, 2., 2.82843, 4., 5.65685, 8., 11.3137, 16., 22.6274, \
32.}

Now we plot these values using ListPlot. Ignore for now the x-tick points, which are just the number of the data point in the List.

In[242]:= plot1 = ListPlot@samplePoints

To better visualize the function, it's often helpful to connect the sample points. We use ListLinePlot. Another way is to use ListPlot's option, Joined -> True). If you look closely, you can see that the sample points are joined by line segments; our plot is not a continuous curve. Again, ignore the abscissa values.

plot2 = ListLinePlot@samplePoints

If we wished, we could combine the plots with Show--one of the most commonly-used plot techniques.

Show[plot1, plot2]

And this is similar to what Plot does--sampling a function and plotting the connected sample points--but of course in a much more sophisticated way. For one thing, Plot samples more frequently where it thinks it might need to, such as where the value of the function changes more frequently, maybe with several tries, and automatically connects the sample points in a continuous curve. Using Plot the x-axis values are correctly labeled.

Plot[f[2, t], {t, -5, 5}]

But notice that Plot is not showing the full domain and range that we asked for. By default it tries to, in effect, show the most important part of a function by excluding what it thinks might be outliers. Two options allow you to override this behavior. PlotRange->Full tells Plot to not clip the y-axis. PlotRange->All tells Plot to include all domain points and their y-values. I usually use PlotRange->All.

Plot[f[2, t], {t, -5, 5}, PlotRange -> All]

If you wish, reveal the sample points with the Option, Mesh->All. When you don't understand what Plot is doing, revealing the Mesh is a good first step. And then as I said above, you can compare Plot's sample points to yours, which might illuminate what's going on. There are more Mesh options to allow you to control how Plot uses its sample points.

Plot[f[2, t], {t, -5, 5}, Mesh -> All, PlotRange -> All]

Next: Using Plot to Plot More than One Function