Posts Tagged ‘monte carlo integration’

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:

I=∫abf(x)dx
I≈∑i=1nf(xi)p(xi)

xi is a random value within the range [a,b]

p(x) represents the distribution of random values, for a uniform distribution:

p(x)=1b–a

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();
}

0