Hello all,

I have the mcintegrate.h algorithm working correctly and I have run the examples given in the NR3 book, however I am wondering how to call the mcintegrate routine for a simple double integral.

Here is a simple double integral example taken from this website.

Example 6

http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html

How would you call mcintegrate to perform the Monte Carlo integration of that simple function?

Any help would be greatly appreciated.

Thanks
Matt

...

Here is a simple double integral example taken from this website.

Example 6

http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html

How would you...

//
// Numerical Recipes mcintegrate to find approximate
// integral of f(x,y) = 4 - x^2 - y^2
// over a rectangular region bounded by [0,5/4][0,5/4]
//
// davekwx7
//
#include "../code/nr3.h"
#include "../code/ran.h"
#include "../code/mcintegrate.h"
//
// This is the function to be integrated:
// z = 4 - x^2 - y^2
//
VecDoub func(const VecDoub & x)
{
VecDoub f(1);
f[0] = 4.0 - x[0]*x[0] - x[1]*x[1];
return f;
}

//
// The region is a simple rectangle, and the
// simple definition of minimum and maximum values
// of x and y assures that the point will be
// in the box.
//
// Therefore: always return true
//
Bool region(const VecDoub & x)
{
return true;
}
int main()
{
Int numsteps;
Int numruns = 3;
// A simple rectangular region:
Doub xmin = 0.0;
Doub ymin = 0.0;
Doub xmax = 5.0/4.0;
Doub ymax = 5.0/4.0;

VecDoub lower(2);
VecDoub upper(2);
lower[0] = xmin;
lower[1] = ymin;
upper[0] = xmax;
upper[1] = ymax;

//
// See the derivation of the analytical result at
//http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html
//
Doub actual = 1775.0 / 384.0;
cout << fixed;
//
// Different random sequence each time you run it
// (Assuming the runs are at least a second apart.
Int ranseed = time(0);

//
// For debugging and comparison: same seed each time
//
//Int ranseed = 87654321;

cout << "Enter the number of Monte Carlo steps: ";
cin >> numsteps;
if ((!cin) || (numsteps < 1)) {
cout << "Invalid entry. Must be a positive integer." << endl;
exit(EXIT_FAILURE);
}

cout << endl
<< "Analytical result = " << actual
<< endl << endl;

cout << "Number of Monte Carlo steps = " << numsteps
<< endl << endl;

for (Int i = 0; i < numruns; i++) {
ranseed += numsteps*i; // Different seed each run
MCintegrate mcint(lower, upper, func, region, NULL, ranseed);
mcint.step(numsteps);
mcint.calcanswers();
cout << "Run " << i+1 << endl
<< " Approximate value = " << mcint.ff[0] << endl
<< " Estimated abs err = " << mcint.fferr[0] << endl
<< " Actual absolute err = " << fabs(mcint.ff[0] - actual)
<< endl << endl;
}

return 0;
}


A few runs:

Enter the number of Monte Carlo steps: 100

Analytical result = 4.622396

Number of Monte Carlo steps = 100

Run 1
Approximate value = 4.657491
Estimated abs err = 0.098498
Actual absolute err = 0.035095

Run 2
Approximate value = 4.628399
Estimated abs err = 0.105627
Actual absolute err = 0.006004

Run 3
Approximate value = 4.419497
Estimated abs err = 0.108110
Actual absolute err = 0.202899



Enter the number of Monte Carlo steps: 1000

Analytical result = 4.622396

Number of Monte Carlo steps = 1000

Run 1
Approximate value = 4.679553
Estimated abs err = 0.031540
Actual absolute err = 0.057157

Run 2
Approximate value = 4.670782
Estimated abs err = 0.032336
Actual absolute err = 0.048387

Run 3
Approximate value = 4.652503
Estimated abs err = 0.033351
Actual absolute err = 0.030107



Enter the number of Monte Carlo steps: 10000

Analytical result = 4.622396

Number of Monte Carlo steps = 10000

Run 1
Approximate value = 4.617056
Estimated abs err = 0.010308
Actual absolute err = 0.005340

Run 2
Approximate value = 4.632975
Estimated abs err = 0.010276
Actual absolute err = 0.010579

Run 3
Approximate value = 4.620658
Estimated abs err = 0.010269
Actual absolute err = 0.001738



Regards,

Dave

Dave,

Thanks for your help. Everything works well.

Matt

Dave,

Would you be able to show how to implement a timing loop for the mcintegate() using the <ctime> library?

Thanks
Matt

... how to implement a timing loop... Let's forget mcintegrate for now (I'll get back to it in a minute or two).

Here's the thing about delay loops in standard C or C++:

Many people seem to think there is a standard library "sleep()" or, maybe, "Sleep()" function, but, in fact, there isn't. There are a couple of Posix "sleepy-time" functions that the good GNU gcc folks make available, but Borland and Microsoft don't have them. On the other hand, compilers that are used in Windows environments usually have a <windows.h> header that gives a prototype for a "Sleep()" function whose argument is in milliseconds. (You can include <windows.h> even if you are writing a command-line program that doesn't use the Windows API.)

GNU compilers on Linux don't have <windows.h> but they have a function called "usleep()" whose prototype is in <unistd.h> Its argument is in microseconds.

So, here is just about the simplest reasonably portable and reasonably efficient "delay loop" program that can be compiled with my GNU compilers on Linux and also with GNU, Microsoft or Borland compilers on my Windows platform. Try it and see if it gives you an output line every second or so. (See Footnote.)


#include <iostream>
#ifdef _WIN32
#include <windows.h>
#else
#ifdef __GNUC__
#include <unistd.h>
// Sleep for milliseconds
void Sleep(unsigned t)
{
usleep(1000*t);
}
#else
#error "Not _WIN32 and not __GNUC__;don't know how to sleep"
#endif
#endif

using namespace std;

int main()

{
for (int i = 0; i < 10; i++) {
cout << i+1 << endl;
Sleep(1000); // 1000 ms = 1 second
}
return 0;
}


You can change the argument for the Sleep function for more or less delay between output lines. If that doesn't work, then stop here and tell me what compiler and operating system you are using.



Now, as far as mcintegrate delay loops:

Are you saying that you want to print out values from the steps of mcintegrate() in some kind of delay loop or what?

For example:

// The following could be put in its own header,
// maybe Sleep.h or some such thing
//
/*---------Start of special Sleep stuff------------*/
#ifdef _WIN32
#include <windows.h>
#else
#ifdef __GNUC__
#include <unistd.h>
// Sleep for milliseconds
void Sleep(unsigned t)
{
usleep(1000*t);
}
#else
#error "Not _WIN32 and not __GNUC__;don't know how to sleep"
#endif
#endif
/*----------End of special Sleep stuff-------------*/

//
// Numerical Recipes mcintegrate to find approximate
// integral of f(x,y) = 4 - x^2 - y^2
// over a rectangular region bounded by [0,5/4][0,5/4]
//
// davekwx7
//
#include "../code/nr3.h"
#include "../code/ran.h"
#include "../code/mcintegrate.h"
//
// This is the function to be integrated:
// z = 4 - x^2 - y^2
//
VecDoub func(const VecDoub & x)
{
VecDoub f(1);
f[0] = 4.0 - x[0]*x[0] - x[1]*x[1];
return f;
}

//
// The region is a simple rectangle, and the
// simple definition of minimum and maximum values
// of x and y assures that the point will be
// in the box.
//
// Therefore: always return true
//
Bool region(const VecDoub & x)
{
return true;
}
int main()
{
Int numsteps;
Int numruns = 3;
// A simple rectangular region:
Doub xmin = 0.0;
Doub ymin = 0.0;
Doub xmax = 5.0/4.0;
Doub ymax = 5.0/4.0;

VecDoub lower(2);
VecDoub upper(2);
lower[0] = xmin;
lower[1] = ymin;
upper[0] = xmax;
upper[1] = ymax;

//
// See the derivation of the analytical result at
//http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html
//
Doub actual = 1775.0 / 384.0;
cout << fixed;
//
// Different random sequence each time you run it
// (Assuming the runs are at least a second apart.
//Int ranseed = time(0);

//
// For debugging and comparison: same seed each time
//
Int ranseed = 87654321;

cout << "Enter the number of Monte Carlo steps: ";
cin >> numsteps;
if ((!cin) || (numsteps < 1)) {
cout << "Invalid entry. Must be a positive integer." << endl;
exit(EXIT_FAILURE);
}

cout << endl
<< "Analytical result = " << actual
<< endl << endl;

cout << "Number of Monte Carlo steps = " << numsteps
<< endl << endl;

MCintegrate mcint(lower, upper, func, region, NULL, ranseed);
cout << " Step Approximate Value" << endl;
for (Int n = 0; n < numsteps; n++) {
mcint.step(1);
mcint.calcanswers();
cout << setw(5) << n+1
//<< ": Approximate value = " << mcint.ff[0] << endl;
<< setw(16) << mcint.ff[0] << endl;
Sleep(1000); // 1000 ms = 1 second
}
cout << endl << "Final step " << numsteps << endl
<< " Analytical result = " << actual << endl
<< " Approximate value = " << mcint.ff[0] << endl
<< " Estimated abs err = " << mcint.fferr[0] << endl
<< " Actual absolute err = " << fabs(mcint.ff[0] - actual)
<< endl << endl;

return 0;
}


Output lines every second or so:

Enter the number of Monte Carlo steps: 10

Analytical result = 4.622396

Number of Monte Carlo steps = 10

Step Approximate Value
1 4.054708
2 4.247127
3 4.495900
4 4.433476
5 4.305183
6 4.618923
7 4.595402
8 4.727343
9 4.700996
10 4.646083

Final step 10
Analytical result = 4.622396
Approximate value = 4.646083
Estimated abs err = 0.225694
Actual absolute err = 0.023687


Or some such thing.


Regards,

Dave

Footnote: If portability is not an issue, and you don't want to mess with that conditionally compiled "special Sleep stuff" (or if it doesn't work with your compiler) try the following:

If your compiler has <windows.h>, just include that header and use Sleep(). (Argument is milliseconds.)

If your non-Windows compiler has <unistd.h>, then you can just include that header and use usleep() directly in your program. (Argument is microseconds.)

Try it and see ... (See Footnote.)

Yes, that works fine.

Are you saying that you want to print out values from the steps of mcintegrate() in some kind of delay loop or what?

What I want to do is to have the program return how much time it required to perform the Monte Carlo integration. I know I can do this to a small degree of accuracy using a watch but I know that the computer can do a much better job especially if I forget to wear my watch that day.

Thanks
Matt

...have the program return how much time it required to perform the Monte Carlo integration....
Well that's not quite the same as implementing a timing loop (actually it's not even close to the same), which is what you asked for.

So: Here we go again. It depends...

Rather than wasting the bandwidth and boring you with needless details, please tell us:

What compiler? What operating system?


Regards,

Dave

Well that's not quite the same as implementing a timing loop (actually it's not even close to the same), which is what you asked for.

Sorry about the mix up.

I am using the following

Operating System - Microsoft Windows XP Professional Version 2002 Service Pack 3

Compiler - Visual C++ 2008 Express Edition. Version 3.5 SP 1

Thanks
Matt

Visual C++ 2008...

For starters you can try the clock() function (which is a standard library function for C and C++).

This gives elapsed time since the process started. Even though it may seem that the resolution is in microseconds, the granularity depends on the operating system so the minimum time difference is on the order of a few tens of milliseconds. Also note that the reading depends on what else is going on in the system at that time (other processes). If you need better resolution, see Footnote.

Anyhow, you can try something like this:


//
// Numerical Recipes mcintegrate to find approximate
// integral of f(x,y) = 4 - x^2 - y^2
// over a rectangular region bounded by [0,5/4][0,5/4]
//
// davekwx7
//
#include "../code/nr3.h"
#include "../code/ran.h"
#include "../code/mcintegrate.h"
//
// This is the function to be integrated:
// z = 4 - x^2 - y^2
//
VecDoub func(const VecDoub & x)
{
VecDoub f(1);
f[0] = 4.0 - x[0]*x[0] - x[1]*x[1];
return f;
}

//
// The region is a simple rectangle, and the
// simple definition of minimum and maximum values
// of x and y assures that the point will be
// in the box.
//
// Therefore: always return true
//
Bool region(const VecDoub & x)
{
return true;
}
int main()
{
Int numsteps;
// A simple rectangular region:
Doub xmin = 0.0;
Doub ymin = 0.0;
Doub xmax = 5.0/4.0;
Doub ymax = 5.0/4.0;

VecDoub lower(2);
VecDoub upper(2);
lower[0] = xmin;
lower[1] = ymin;
upper[0] = xmax;
upper[1] = ymax;

//
// See the derivation of the analytical result at
//http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html
//
Doub actual = 1775.0 / 384.0;
//
// Different random sequence each time you run it
// (Assuming the runs are at least a second apart.
Int ranseed = time(0);

//
// For debugging and comparison: same seed each time
//
//Int ranseed = 87654321;

cout << "Enter the number of Monte Carlo steps: ";
cin >> numsteps;
if ((!cin) || (numsteps < 1)) {
cout << "Invalid entry. Must be a positive integer." << endl;
exit(EXIT_FAILURE);
}

cout << endl
<< "Analytical result = " << actual
<< endl << endl;

cout << "Number of Monte Carlo steps = " << numsteps
<< endl << endl;

MCintegrate mcint(lower, upper, func, region, NULL, ranseed);

clock_t start_time = clock();
mcint.step(numsteps);
mcint.calcanswers();
clock_t end_time = clock();

double elapsed_time = end_time - start_time;
cout << "Elapsed time = "
<< elapsed_time / CLOCKS_PER_SEC
<< " seconds." << endl << endl;

cout << fixed;
cout << " Approximate value = " << mcint.ff[0] << endl
<< " Estimated abs err = " << mcint.fferr[0] << endl
<< " Actual absolute err = " << fabs(mcint.ff[0] - actual)
<< endl << endl;

return 0;
}


Output for a run on my creaky old Windows XP platform:

Enter the number of Monte Carlo steps: 1000000

Analytical result = 4.6224

Number of Monte Carlo steps = 1000000

Elapsed time = 0.922 seconds.

Approximate value = 4.622233
Estimated abs err = 0.001030
Actual absolute err = 0.000163


Regards,

Dave

Footnote: Most Windows compiler users can get microsecond-resolution timing information with the Windows API functions QueryPerformanceCounter() and QueryPerformanceFrequency(). I can give an example if you are interested.

Output can look like

Enter the number of Monte Carlo steps: 100

Analytical result = 4.6224

Number of Monte Carlo steps = 100

Elapsed time = 0.000100851 seconds

Approximate value = 4.688166
Estimated abs err = 0.103545
Actual absolute err = 0.065770

Dave, that works perfect.

I can give an example if you are interested.

Sure, if you have the time, I would be interested.

Thanks
Matt

...interested...

//
// Numerical Recipes mcintegrate to find approximate
// integral of f(x,y) = 4 - x^2 - y^2
// over a rectangular region bounded by [0,5/4][0,5/4]
//
// davekwx7
//
// Recent Microsoft, Borland, GNU Windows compilers can
// (probably) use high-resolution timer if you include
// <windows.h>
//
#include <windows.h>

#include "../code/nr3.h"
#include "../code/ran.h"
#include "../code/mcintegrate.h"

//
// This is the function to be integrated:
// z = 4 - x^2 - y^2
//
VecDoub func(const VecDoub & x)
{
VecDoub f(1);
f[0] = 4.0 - x[0]*x[0] - x[1]*x[1];
return f;
}

//
// The region is a simple rectangle, and the
// simple definition of minimum and maximum values
// of x and y assures that the point will be
// in the box.
//
// Therefore: always return true
//
Bool region(const VecDoub & x)
{
return true;
}
int main()
{
LARGE_INTEGER start_pc;
LARGE_INTEGER end_pc;
LARGE_INTEGER freq_pc;

Int numsteps;

// A simple rectangular region:
Doub xmin = 0.0;
Doub ymin = 0.0;
Doub xmax = 5.0/4.0;
Doub ymax = 5.0/4.0;

VecDoub lower(2);
VecDoub upper(2);
lower[0] = xmin;
lower[1] = ymin;
upper[0] = xmax;
upper[1] = ymax;

//
// See the derivation of the analytical result at
//http://math.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_7.html
//
Doub actual = 1775.0 / 384.0;
//
// Different random sequence each time you run it
// (Assuming the runs are at least a second apart.
Int ranseed = time(0);

//
// For debugging and comparison: same seed each time
//
//Int ranseed = 87654321;


//
// First make sure your system supports the high resolution
// timer functions.
//
// You could use FormatMessage with GetLastError to
// print out the specific cause of failure, but
// I'll just bail out with a minimal message.
//
if (!QueryPerformanceFrequency(&freq_pc)) {
cerr << "QueryPerformanceFrequency() failed." << endl;
exit(EXIT_FAILURE);
}
if (!QueryPerformanceCounter(&start_pc)) {
cerr << "QueryPerformanceCounter() failed." << endl;
exit(EXIT_FAILURE);
}

cout << "Enter the number of Monte Carlo steps: ";
cin >> numsteps;
if ((!cin) || (numsteps < 1)) {
cout << "Invalid entry. Must be a positive integer." << endl;
exit(EXIT_FAILURE);
}

cout << endl
<< "Analytical result = " << actual
<< endl << endl;

cout << "Number of Monte Carlo steps = " << numsteps
<< endl << endl;

MCintegrate mcint(lower, upper, func, region, NULL, ranseed);

QueryPerformanceCounter(&start_pc);
mcint.step(numsteps);
mcint.calcanswers();
QueryPerformanceCounter(&end_pc);

double diff_ct = end_pc.QuadPart - start_pc.QuadPart;

cout << "Elapsed time = " << diff_ct / freq_pc.QuadPart << " seconds."
<< endl << endl;

cout << fixed;
cout << " Approximate value = " << mcint.ff[0] << endl
<< " Estimated abs err = " << mcint.fferr[0] << endl
<< " Actual absolute err = " << fabs(mcint.ff[0] - actual)
<< endl << endl;

return 0;
}


Output for a thousand steps takes about a millisecond on my very modest Windows XP workstation:

Enter the number of Monte Carlo steps: 1000

Analytical result = 4.6224

Number of Monte Carlo steps = 1000

Elapsed time = 0.000954311 seconds.

Approximate value = 4.660326
Estimated abs err = 0.032790
Actual absolute err = 0.037930

Thanks Dave. I always love examples.

Matt