5

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?

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
johny
  • 355
  • 3
  • 7

4 Answers4

4

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}
Community
  • 1
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
2

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

Pankaj Sejwal
  • 2,063
  • 14
  • 23
2

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}
george2079
  • 38,913
  • 1
  • 43
  • 110
0

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.

Pankaj Sejwal
  • 2,063
  • 14
  • 23
  • 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. – Pankaj Sejwal 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. – Pankaj Sejwal Nov 15 '13 at 20:05