Trivial example: Monte Carlo integration of the normal distribution

This is a very simple (and not in itself useful) example of using floating point random numbers in practice. To goal is to evalute the integral of the standard normal distribution:

The approach taken is the so-called Monte-Carlo integration where randomly drawn numbers are used to explore the function to be integrated.

The Monte-Carlo algorithm

The algorithms is based on drawing numbers according to the standard normal distribution and simply working out the fraction of those that are within the limits of the integration:

This approach is quite different from the alternative of drawing uniformly distributed random numbers and then weigthing each draw with the value of the function, in that we are automatically going to preferentially sample the most important part of the distribution.

The program

The implementation of the above algorithm is very straightforward using the boost normal random number generator already desribed (Normally distributed random numbers):

#include <iostream>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>


/** Integrates the standard normal distribution between xlow and xhigh
    
    \param n is the number of samples (or draws) to make
 */
double mcinteg(double xlow,
               double xhigh,
               size_t n)
{
  boost::mt19937 igen;
  boost::variate_generator<boost::mt19937, boost::normal_distribution<> >
    gen(igen, 
        boost::normal_distribution<>());

  size_t k=0;
  for(size_t i=0; i<n; ++i)
  {
    const double x=gen();
    if (x>= xlow and x<xhigh)
      ++k;
  }
  return (double)k/n;
}

int main(void)
{
  const double xlow=-1.0;
  const double xhigh=1.0;

  std::cout<<"MC Result: "
           <<mcinteg(xlow, xhigh, 1000000)
           <<std::endl
           <<"Expected: "
           <<0.5*(erf(xhigh/sqrt(2.0))-erf(xlow/sqrt(2.0)))
           <<std::endl;
}

The core of the program is the for loop in mcinteg which draws n random numbers, checks if each one is within the limits of integration, and if it is, adds one to the running total of draws that were within the limits.

The other signficant of the program is the line which computes the expected value of the integral using the library implementation of the error function.

The output of the program is:

MC Result: 0.682773
Expected: 0.682689

And as you can see, with samples the computed value is good to about 4 decimal places.