[Exercise] Approximating Pi - Monte Carlo Analysis

Ra8

this exrcise was inspired from the one rahmu posted
I'm taking introduction to programming at university, our professor gave us a similar exercice (approximate pi using monte carlo analysis)

this is the exercise:

Monte Carlo simulations are techniques that consist of choosing sample experiments at random from a large set and then making deductions on the basis of the probabilities estimated from the result of these experiments. In this problem, you are to implement a Monte Carlo method for estimating the value of π. Consider the first quadrant of a unit circle centered at (0,0). This quarter circle lies inside a unit square. A point with coordinates (x,y) is inside the quarter circle iff x2 + y2 <= 1. The area of the quarter circle region can be estimated by picking, at random, points (x,y) that lie in the unit square, and for each point determining whether the point lies in the region. The fraction of points that fall in the region should give an estimate of π/4 (ratio of the area of the region and the area of the enclosing unit square). Multiplying by 4 gives an estimate of π. Write a program that takes a command line integer parameter N and prints an estimate of π using N random points as described above.

my solution in C++:

#include <iostream>
#include <stdlib.h>
#include <ctime>
using namespace std;

int main()
{
    srand(time(NULL));
    float rand_x, rand_y;
    int N, pointsRegion=0;
    cout << "Enter number of points: ";
    cin >> N;

    for(int i=0; i<N; i++){
        rand_x= rand() * 1.0 / RAND_MAX;
        rand_y= rand() * 1.0 / RAND_MAX;
        if((rand_x*rand_x+rand_y*rand_y)<=1)
        {
            pointsRegion++;
        }
    }
    cout << "\n\nAn estimatimation of pi with " << N << " random points is: " << (pointsRegion*4.0/N) <<"\n\n";
    return 0;
}

geek

With[{n = 10^6}, 4 Length@Select[RandomReal[{0, 1}, {n, 2}], Norm@# <= 1.0 &] / n // N]

Joe

I wanted to do the exercise in C, because it's been a while since I used the language. I ended up with something extremely similar to what Ra8 did. So I decided to add the possiblity to run the test multiple times and averaging the results.

In retrospect that was generally useless because I could've just increased the value of N; but all in all it got me manipulating some C so good experience.

#include <stdlib.h>
#include <stdio.h>
#include <time.h>


int main (int argc, char **argv)
{
    int i,j; /* Loop counter */
    double x,y; /* Tmp coordinates */
    int hit=0; /* Positive hits */
    double s=0; /* sum for later averaging */
        
    double *run = NULL; /* A list that holds the result of each run */

    int n=1000000; /* Number of iteration per run */
    int runnbr=10; /* Number of runs */


    /* Initializing the srand module, and the run dynamic list */
    srand(time(NULL));
    run = (double *)malloc (sizeof(double) * runnbr);


    /* The run loop */
    for (j=0; j<runnbr; ++j) {

        hit=0; /* Initializes the hit counter before each run */

        /* The following is repeated n times in each run */
        for (i=0; i<n; ++i) {
            x=rand()/(double)RAND_MAX;
            y=rand()/(double)RAND_MAX;

            if ((x*x) + (y*y)<=1)
                ++hit;
        }

        /* Store the result of each run in the run array */
        run[j]=(double)hit/n*4.0;
    }

    /* Averaging the result of all the runs */
    for (j=0; j<runnbr; ++j)
        s += run[j];
            
    printf("%f\n", s/runnbr);

    
    return EXIT_SUCCESS;
}

Joe

Also, the code in Scheme (you can guess I'm bored out of my mind at the office today :P )

(define square
  (lambda (x) (* x x)))
 
(define (f-iter n count hit)
  (cond ((> count n)(/ (* hit 4.0) n))
        ((<= (+ (square (random 1.0))
                (square (random 1.0)))
             1) (f-iter n (+ count 1) (+ hit 1)))
        (else (f-iter n (+ count 1) hit))))
 
 
(define mc-pi
  (lambda (n)
    (f-iter n 0 0)))

Result:

> (display (mc-pi 1000000))
3.141584