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

]