Help with monte carlo integration in C++

[assassin_696]

I'm trying to write a program that calculates the value of ln(2) using monte carlo methods in C++. Obviously ln(2) can be calculated from the integral of 1/x between the limits 2 and 1, but I'm having some trouble implementing that into code.

 

So far I've got:

 

//Program uses Monte Carlo (random number generator) methods
//to calculate a value for ln(2)

#include <iostream>
#include <cstdlib>
#include <cstdio>

using namespace std;

#define ranf() ??????????

int main(int argc, char* argv[])
{
   int outcome, N=10, count_in=0, seed=30;
   double fraction_in;

   if(argc>1)   {
   sscanf( argv[1], "%d", &N ) ; // puts first command line argument in N
   }
   if(argc>2) {
   sscanf( argv[2], "%d", &seed ); // puts the 2nd argument in seed
   }
   cout << " # " << argv[0] << " N= " << N << " seed= " << seed << endl;

   srand(seed); // initialises the random number generator

   for(int n=1; n<=N; n++) {
   double x =ranf();    // ranf() needs to return a real number in [0,1]
   double y =ranf();
   outcome = (x*y < 1.0 ) ? 1 :0;
   count_in += outcome;
   fraction_in = static_cast<double>(count_in)/n; //casts integer variables for correct division
   cout << "Location" << outcome << "\t" << x << "\t" << y << "\t"
   << count_in << "\t" << n << "\t" << fraction_in <<endl;
   }
   return 0;
}

 

The trouble I'm having is putting the limits of the integration in, which I think is done in the ranf() function where I put a load of question marks. Looking at a similar example for calculating the value of pi by testing whether or not a point is in the unit circle, ranf is

 

#define ranf() \
((double)rand()/(1.0+(double)0x7fffffff))

 

Which I think is where the limits information is obtained, but I'm not sure. I also want to make sure that my outcome expression, i.e.

 

 outcome = (x*y < 1.0 ) ? 1 :0; 

 

Is correct, since that's where the integral equation information is contained, obviously if I screwed that up it'd be calculating the wrong integral.

 

Any help would be appreciated, I'm pretty new to all this so not sure how easy that code is to read for more experienced people but obviously I can clarify anything.

[Markup]

Updating, ignore

 

http://www.cplusplus.com/reference/std/complex/log/

http://www.cplusplus.com/reference/clibrary/cstdlib/rand/

 

The second link being the most relevant. It shows you how to setup a seeded random generator. It also explains how to limit the range using modulo.

 

replaced!

 

Enjoy

 

Notes:

 

Observe indentation, although opinion based ;p

The declaration of x and y are outside of the for loop to stop performance impact of declaring the variable on each iteration of the loop.

 

outcome = (x*y < 1.0 ) ? 1 : 0;

 

If condition is true, return the first value(left side of : ), if false return the second value(right side of : )

Edited February 19, 2011 by Markup

[Veiva]

In order to generate a number in the range of [0, 1], the usual best method is to divide the return value of rand by the value of RAND_MAX. Like so:

 

#define randf() (rand() / (double)RAND_MAX)

 

That will assuredly give you a value between 0 and 1. You can verify the math yourself using a set using an example implementation of rand that has a maximum return value of 255. For any value between 0 and 255, divide it by 255. You'll get 0, 0.004, 0.008, ..., 1.

[dsavi]

Wow I wish I understood this.

[Markup]

Wow I wish I understood this.

 

I understand the monte carlo concept, just unsure on how to implement the finer details. In layman's terms,

You place a box around the area of a graph you want to integrate, then using a random sampling method, you test points in the box. You keep track of how many points tested for and how many are under the line. Using this data you can approximate areas or expressions.

 

eg/

An illustration of Monte Carlo integration. In this example, the domain D is the inner circle and the domain E is the square. Because the square's area can be easily calculated, the area of the circle can be estimated by the ratio (0.8) of the points inside the circle (40) to the total number of points (50), yielding an approximation for pi/4 approx 0.8

 

(PI*(d^2))/4 = (40/50)*(d^2)

 

=> PI/4 ~ (40/50)

 

The calculated approx area being, (40/50)*(d^2) = 3.2

[Markup]

Think I got it

 

 

//Program uses Monte Carlo (random number generator) methods
//to calculate a value for ln(2)

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

#define myRand() ((double)rand() / (double)RAND_MAX);

int main(int argc, char* argv[])
{
int seed = (int)time(NULL);
int n=100000;
double x;
double y;

int r = 0;
double ratio = 0.0;

srand(seed);

for(int i = 0; i < n; i++)
{
	x = myRand();
	y = myRand();

	x += 1.0;

	r += (y*x <= 1.0) ? 1 : 0;

	//uncomment below to enter the matrix!
	//ratio = (double)r/(double)i;
	//std::cout << "Iteration: " << i << "\tArea: " << ratio*2.0 << "\tRatio: " << ratio << "\tx: " << x << "\ty: " << y << std::endl;
}
ratio = (double)r/(double)i;
std::cout << "Ratio: " << ratio << "\tArea: "<< ratio*2.0 << std::endl;

return 0;
}