Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

This question already has an answer here:

I wrote the following code to find the area of the circle $x^2+y^2=1$ using monte carlo simulation :

McArea[Num_] := Module[{hit, miss, index, x, y}, 
  hit = 0; miss = 0;
  For[index = 1, index <= Num, index = index + 1,
    x = Random[Real, {-1, 1}];
    y = Random[Real, {-1, 1}];
    If[y <= Sqrt[1 - x^2], hit = hit + 1, miss = miss + 1];];
  Return[(hit/Num) 4];]

McArea[100] is around 3.56, and the larger I make Num, the result gets larger, but only a bit. It doesn't seem to be the correct because my output is too far from 3.14.

Is there something wrong with the code?

share|improve this question

marked as duplicate by Szabolcs, m_goldberg, ybeltukov, Artes, rm -rf Nov 15 '13 at 23:37

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

    
Are you using Mathematica V6? –  Kuba Nov 15 '13 at 11:18
    
@Kuba No Mathematica 8.0. –  johny Nov 15 '13 at 11:18
4  
I would use If[Abs[y] <= ... or two Ifs, one for positive y and one for negative y. –  Peltio Nov 15 '13 at 11:21
    
Then You can use RandomReal. Also take a look at loops in Mathematica what may help you in improving the performance of the code. p.s. why are you counting miss? –  Kuba Nov 15 '13 at 11:27
2  
Look at RegionPlot[y <= Sqrt[1 - x^2], {x, -1, 1}, {y, -1, 1}] –  Simon Woods Nov 15 '13 at 19:12

4 Answers 4

up vote 3 down vote accepted

Your simulation's convergence to an incorrect result is because your comparison accepts any negative y, thus inflating the result. You are in essence computing the area of the union of the semi-unit circle that lies above the x-axis and the rectangle with top-left corner at {-1,0} and lower-right corner at {1,-1}, which is 2 + π/2 = 3.5708.

This is the minimal repair to your code.

McArea[Num_] := 
  Module[{hit, miss, index, x, y}, 
    hit = 0; miss = 0;
    For[index = 1, index <= Num, index = index + 1, 
      x = Random[Real, {-1, 1}];
      y = Random[Real, {-1, 1}];
      If[Abs[y] <= Sqrt[1 - x^2], hit = hit + 1, miss = miss + 1];];
    Return[(hit/Num) 4];]
Timing[N @ McArea[100000]]
 {0.590944, 3.13596}

Now let's explore some improvements. The first still uses a For-loop, but with a slight better comparison expression and some improved Mathematica practice.

mcArea2[num_] :=
  Module[{hit, index, x, y},
    For[hit = 0; iindex = 1, index <= num, ++index, 
      {x, y} = RandomReal[{-1, 1}, 2];
      If[x^2 + y^2 <= 1, ++hit]];
    4. hit/num]
Timing[mcArea2[100000]]
{0.571963, 3.14488}

Note the use of lowercase initial letters on all user defined identifiers. Using uppercase initial letters is bad practice because it can get you into conflicts with Mathematica's predefined, built-in variables. Also, note the use ++ and the removal of the unnecessary Return.

Do is better than For in this case since we don't need any of fine control over iteration that For provides.

mcArea3[num_] :=
  Module[{hit, x, y},
    hit = 0;
    Do[{x, y} = RandomReal[{-1, 1}, 2]; If[x^2 + y^2 <= 1, ++hit], {num}];
    4. hit/num]
Timing[mcArea3[100000]]
{0.476420, 3.14056}

Edit

People have started to submit functional programming answers where the connection to the original question is no longer obvious, so I will add mine.

mcArea4[num_] :=
  4. Plus @@ Table[Boole[Norm[RandomReal[{-1, 1}, 2]] <= 1], {num}]/num
Timing[mcArea4[100000]]
{0.032864, 3.14664}
share|improve this answer

This one will give correct result,

   McArea[Num_] := Module[{hit, miss, index, x, y, ar}, hit = 0; miss = 0;
      For[index = 1, index <= Num, index = index + 1,
       x = RandomReal[{-1, 1}];
       y = RandomReal[{-1, 1}];
       ar = x^2 + y^2;
       If[ar <= 1, hit = hit + 1, miss = miss + 1]];
      Return[(hit/Num) 4];]

McArea[1000000] // N

3.14101

share|improve this answer

better still..

      4 Length@
        Select[Table[RandomReal[{-1, 1}, 2], {100000}], Norm[#] < 1 &]/
        100000 // N // Timing

      {0.171601, 3.14396}

or

      4 Count[ Table[RandomReal[{-1, 1}, 2], {100000}] , x_ /; Norm[x] <1]/
          100000 // N // Timing

      {0.234002, 3.14192}
share|improve this answer

Just a minor change in your range does it,

McArea[Num_] := Module[{hit, miss, index, x, y}, hit = 0; miss = 0;
  For[index = 1, index <= Num, index = index + 1, 
   x = RandomReal[{-1.31, 1.31}];
   y = RandomReal[{-1.31, 1.31}];
   If[y <= Abs[Sqrt[1 - x^2]], hit = hit + 1, miss = miss + 1];];
  Return[(hit/Num) 4];]

3.11298

You can hit and try more to get more appropriate results.

share|improve this answer
1  
Somehow, I do not think that this was the problem. The problem was that the simulation was considering all points under the upper semi-circonference, not excluding those 'outside' the lower semi-circonference. That is why it was returning value in excess of the actual solution. Using random numbers bigger than 1, as you do, makes the square root complex, thus requiring an Abs for the ordering, but what is the use of that? –  Peltio Nov 15 '13 at 18:38
    
@Peltio In my opinion it is helping to take in more points that satisfy the constraint of being less than 1 though they are not.Their Abs[] will make them less than one. It's like considering adjacent points outside circle to decide area of circle. –  Rorschach Nov 15 '13 at 19:25
    
In addition to what @Peltio wrote, you might want to do a ListPlot of the "hit" points. it will indicate a picture that is quite likely not what is expected. Here is code that compares the two approaches to this computation. McArea[num_] := Module[ {hit, miss, index, pts, pts1, pts2}, hit = 0; miss = 0; pts = RandomReal[{-1.31, 1.31}, {num, 2}]; pts1 = Select[pts, #[[2]] <= Abs[Sqrt[1 - #[[1]]^2]] &]; pts2 = Select[pts, Norm[#] <= 1.31 &]; lp = ListPlot[pts1, AspectRatio -> 1]; {Length[pts1], Length[pts2]}/num*4.]. Now do McArea[10000] and then check lp. –  Daniel Lichtblau Nov 15 '13 at 19:48
    
@DanielLichtblau : Thanks for this nice depiction, my assumption turned out to be wrong largely but I think this variation will give approx results at least.It is definitely not what OP asked but a different way in my opinion. –  Rorschach Nov 15 '13 at 20:05

Not the answer you're looking for? Browse other questions tagged or ask your own question.