2

Generally I understand the idea of the Monte Carlo method.

However, when I read articles about it, there is always shown an example of calculating pi using a square, into which we insert 1/4th of a circle. Then we start putting randomly points into the square, and because the area of a 1/4 of the circle is pi/4, and the area of the square is exactly 1, we get that pi/4 = number of points in the circle / number of points in the square.

Monte Carlo example from Wikipedia

This simple example had been obvious for me, I'd liked it much, but one day I realized, that this gives correct (with some error of course) results because we know exactly what we want to get. So in this example we know that this square should be 1-sized.

Let us consider we know nothing about the function (in this case the circle) so we don't know which borders to choose. If we mistakenly take a 2x2 square, we won't get good result.

Ok, so we should first find the square. To find a square we need the function values, at least max and min. So I need to search them. But searching increases complexity. If I check each value for each x why just don't perform numerical integration in the same loop?

I understand the example of searching pi is trivial, however for me this does not make sense as we know how to fit data to get desired result.

(Please note I'm not a mathematician, but an engineer so I have some basis in maths, but I'm not good in advanced problems)

Voitcus
  • 193

2 Answers2

1

This yields a binomial distribution for the points inside the circle. The probability that a point is inside the circle is $p = \left(\pi/4\right)/\left(1\times 1\right) = \pi/4$. Then, the average number of points inside the circle is $\overline{n} = Np = N\pi/4$ where $N$ is the total number of points. That means $$ \pi = 4\,{\overline{n} \over N} $$ In a computer, you generate $N$ points inside the square $\left(~N \gg 1~\right)$ and counts how many are inside de circle ( that is $\overline{n}$ ). This little script ( in $\tt C++$ ) generates a Montecarlo $\pi$:

// piForSE0.cc
#include <cstdlib>
#include <iostream>
using namespace std;
typedef unsigned long long ullong;
const ullong N=1000000ULL; // One million
const double RANDMAX1=double(RAND_MAX) + 1.0;

inline double rand01()
{
 return rand()/RANDMAX1;
}

int main()
{
 double x,y;
 ullong pointsInCircle=0;

 for ( ullong n=0 ; n<N ; ++n ) {
     x=rand01();
     y=rand01();
     if ( (x*x + y*y)<1.0 ) ++pointsInCircle; 
 }

 cout<<"PI = "<<(4.0*(double(pointsInCircle)/N))<<endl;
 return 0;
}
Felix Marin
  • 89,464
  • Well, I generally understand this example and how it is performed. I don't understand how to use this method for integration when we don't know which square to use. If instead of 1x1 square we take 2x1 rectangle, this method fails if we don't update the last formula too. – Voitcus Oct 05 '13 at 06:41
  • @Voitcus If you take $2X1$, then $p = (\pi/4)/(2\times 1) = \pi/8$. Then $\overline{n} = Np = N\times(\pi/8)$ and $\pi = 8\times(\overline{n}/N)$. In the $\tt C++$ script, it changes the $y$ to $\tt y = 2.0*rand01()$ since now you generate pairs $(x,y)$ with $0 \leq x < 1$ and $\large 0 \leq y < 2$. – Felix Marin Oct 05 '13 at 06:49
  • I know how and what to change here. But there might be more complex example when I take too large boundary without knowing I should also change the result. Here I take a square that I know it fits exactly the circle. But if I don't know what square to choose, I can't use this method because I don't know what factor use to find out the pi (or integral, because pi/4 is here the area of the 1/4th of the circle) – Voitcus Oct 05 '13 at 07:21
  • OK, I got it. If I calculate integral as ratio of points_in / (points_in + points_out) I need also to multiply this by area of 2x1 rectangle. Here it was 1x1 so the area was 1 and the 1 was dropped and the ratio was exactly pi/4 with or without this 1. Thank you for help. – Voitcus Oct 05 '13 at 07:36
0

I got it.

The problem is that pi/4 value was in fact the integral of the function (in <0, 1>) y = sqrt(1-x^2). And I assumed the value of pi/4 was

pi/4 = points_in / (points_in + points_out)

But this is not true in any case, the proper formula should be

pi/4 = Area_of_rectangle * points_in / (points_in + points_out)

In this case Area_of_rectangle was exactly 1 so it was dropped. If I use other rectangle, eg. 2x1, this needs to be taken into account.

If I know exactly what boundaries to take (ie. how large the rectangle should be to cover the whole function, but this can be any size), I can "calculate" with Monte Carlo the integral of any function.

Voitcus
  • 193