### The Monte Carlo algorithm/method

09Jan07

The Monte Carlo method is a way of solving problems using statistical methods; stochastic technique (which means using random numbers) and probability.

Given some probability, P, that an event will occur in certain conditions, a computer can be used to generate those conditions repeatedly. The number of times the event occurs divided by the number of times the conditions are generated should be approximately equal to P.

Scope: The Monte Carlo methods are applied in various fields of mathematics, physics, economics and etc. Monte Carlo simulation methods are especially useful in studying systems with a large number of coupled degrees of freedom, such as liquids, disordered materials, strongly coupled solids, and cellular structures. The Monte Carlo method is often useful for solving problems in physics and mathematics which cannot be solved by analytical means. They are useful for modeling phenomena with significant uncertainty in inputs, such as the calculation of risk in business.

Calculation of Pi: Let’s make a circle inside a square and let’s take a quadrant out of it. If we take a dice and throw it randomly inside the square, the probability that the dice will end up in the circle area is circle_area/square_area. Now if we replace it with formulas, just for one quadrant, we get: Pi=4*dice_hitting_circle_area/dice_hitting_square_area. So by Pythagoras theorem we just have to take random values for x and y coordinates, square them, sum them up, take the square root out of the sum and if the result is less or equal to 1 then the dice is inside the circle. By repeating this few times we can get the approximation of the number Pi.  #include <iostream>
#include <ctime>
#include <cmath>
using namespace std;

int main()
{
double x, y;            //define x and y coordiantes
int count_inside = 0, till=1000;    //how many tries you want
srand(time(NULL));           //define random values
for(int i=0; i<till; i++)
{
x=pow(double(rand()%1000)/1000,2);    //square random x, x^2
y=pow(double(rand()%1000)/1000,2);    //square random y, y^2
if(sqrt(x+y)<=1)                //check if square root is equal or less to 1
count_inside++;                //if it is then add 1 to count_inside
}
cout << 4*double(count_inside)/till << endl;    //calculate the pi by the formula given above
}
out << 4*double(count_inside)/till << endl; //calculate the pi by the formula given above
}

#### 5 Responses to “The Monte Carlo algorithm/method”

1. 1 Mike

good job!

2. 2 Elton Minetto

Hi, i implemented the Monte Carlo calculation of pi using Python, sockets and threads. It’s in portuguese, but the code is universal
http://www.eltonminetto.net/metodo-de-monte-carlo-distribuido.htm

3. 3 Simon

The real pb is to succeed in doing some random function cause in bigger problems than solving pi, finding a suitable random code is hard.
What is the precision of the pi found ?
I’m doing it in matlab and in fortran with my own sequence to radomize numbers, and i cant’ find nearer than 3.14XXXX then the others (X) are false.
NB: I do the MC with 4000 numbers.

4. 4 Klaus

The precision depends at first on the number of samples and next on the precision of the calculation and the RNG (random number generator). This assumes that the RNG is really random enough, of course.

With 4000 samples, two decimal digits is what you can expect. Try 40.000 and 400.000 samples as see what happens.

As the result is a fraction with the number of samples as denominator, the result is quantized. In your case, the pi/4 intermediate result is quntized to 1/4000, and the pi result is quantized to 1/1000. In other words, it can only provide 3 digits above the noise level.

It may also be noteworthy that the approach can also be tried with a uniform non-random distribution of samples. For example, one can use two nested loops from 0.0 to 1.0 with an increament of 0.01 to generate 10.201 samples (not 10.000, because it’s 101 steps from 0 to 1 here). This approach may be of use when no sufficienty good RNG is available.

5. 5 YT

Just a quick comment about efficiency in the code: There’s no need to take the square root. As long as r^2

*
To prove you're a person (not a spam script), type the security word shown in the picture. ## Social Bookmarking

• Post to Del.icio.us
• Seed This Post
• Send to Digg
• Save on Furl
• Add to My Web 2.0
• Ma.gnolia