Monte Carlo integration
Mar 27 2014 · Math
I was reading a bit about random numbers and remembered that I wrote a simple Monte Carlo integrator in C++ a few years ago. I took a few minutes to cleanup and comment the code, which is presented below.
Monte Carlo integration is simple, but surprisingly powerful:
xi is a random value within the range [a,b]
p(x) represents the distribution of random values, for a uniform distribution:
This presentation by Fabio Pellacini provides a lot more details.
The test code in the main() method computes the integral of sin2(x) in the interval [3,5].
#include <iostream>
#include <cmath>
#include <ctime>
using namespace std;
// Functor base class for encapsulating 1-dimensional function to be integrated
class Function1d
{
public:
virtual double operator()(double x)
{
return x;
}
};
// Functor for sine squared function
class SineSquared : public Function1d
{
public:
double operator()(double x)
{
return pow(sin(x), 2);
}
};
// Monte Carlo integrator class declaration
class MonteCarloIntegrator
{
public:
double Run(int numSamples, Function1d& func, double intervalMin, double intervalMax);
};
// Monte Carlo integrator implementation
// ::Run() method implementation
double MonteCarloIntegrator::Run(int numSamples, Function1d& func, double intervalMin, double intervalMax)
{
double sum = 0;
double div = 1.0 / (double)RAND_MAX;
double intervalScale = 1.0 / (intervalMax-intervalMin);
for(int i=0; i<numSamples; i++)
{
double rnd1 = intervalMin + ( ((double)rand() * div) * (intervalMax-intervalMin) );
sum += (func)(rnd1) / intervalScale;
}
return (1.0/(double)numSamples) * sum;
}
// main() function with test code
void main()
{
MonteCarloIntegrator integrator;
double output = integrator.Run(5000000, SineSquared(), 3, 5);
srand(time(NULL));
std::cout << output << endl;
getchar();
}
