Finding the area of a circle using monte carlo simulation [duplicate]

Question

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?

Answer 1

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}

Answer 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

Answer 3

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}

Answer 4

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.