Monte Carlo Integration

Here we will also introduce Monte Carlo integration.

Recall from basic calculus the mean value theorem for integrals: For $f(x)$ continuous over $[a,b]$, there exists a number, $c$, with $a < c < b$, such that

(1)

\begin{align} {1\over{b-a}} \int _a^b f(x)\,{\mathrm d}x = f(c) \end{align}

Another way to think of this is that the area under the curve is the base width $(b-a)$ times the "average height" $f(c)$, in other words, we can rearrange the above formula into

(2)

\begin{align} \int _a^b f(x)\,{\mathrm d}x = (b-a) \times f(c) \end{align}

---

In some sense, the trapezoidal integration we studied in the previous section is nothing other than a set of midpoint integrations. (Think about it for a minute). Another approach we can take is to try to statistically estimate the midpoint by choosing points at random. In one dimension, this is actually kind of slow and innaccurate. But in higher dimensions it becomes one of the few practical methods of determining integrals.

To state the algorithm concisely, we will:

1. Pick $n$ randomly distributed points $x_1, x_2, x_3, ...x_n$ in the interval of $[a,b]$.
2. Determine the average value of the function $\langle f \rangle = {\frac{1}{n}} \sum_{i=1}^{n} f (x_i)$
3. Compute the approximation to the integral using the mean value theorem $\int _a^b f(x)\,{\mathrm d}x \approx (b-a) \times \langle f \rangle$
4. As a consequence of the Law of Large Numbers and/or the Central limit theorem, we can also estimate the error inherent in this approach, as $\mathrm{Error} \approx (b-a) {\sqrt{ {{ \langle f^2 \rangle - {\langle f \rangle }^2}\over {n}}}$, where $\langle f^2 \rangle = {\frac{1}{n}} \sum_{i=1}^{n} f^2 (x_i)$.

---

Let's consider again the case we studied in the previous section,

(3)

\begin{align} I=\int_0^3 {{x^4 e^x} \over{ (e^x - 1)^2 }} {\mathrm d}x \end{align}

we will use the same function for the integrand as before.

program mc_demo
  implicit none
  integer, parameter :: n=10000

  call random_seed()  ! needed to initialize the random number generator used in MC eval
  call MC_integration(n,3.0)

  contains

    subroutine MC_integration(n,end_val)
      implicit none
      integer, intent :: n
      real, intent :: end_val
      real :: x, integral, integral_err
      real (kind=8) :: f, f2 
      integer :: i

      integral = 0.0

      f= 0.0d0
      f2 = 0.0d0

      do i=1,n
         call random_number(x)
         x=x*end_val  ! random_number only returns uniformly distributed from [0.0, 1.0]

         f = f+integrand(x)
         f2 = f2+ (integrand(x)**2)
      end do

      f=f/n
      f2=f2/n

      integral = (end_val-0.0)*f
      integral_err =  (end_val-0.0)*SQRT((f2 - f**2)/n)
      write (*,*) "# MC integration = ",integral,"+/-",integral_err

    end subroutine MC_integration

    function integrand(x) result (value)
      implicit none
      real :: x
      real :: value

      if (x .lt. 0.00001) then
         x = 0.00001
      end if

      value = (x**4)*EXP(X)/((EXP(X)-1.0)**2)
    end function integrand

  end program mc_demo

---

Let's take a look at the new parts:

integer, parameter :: n=10000

call random_seed()  ! needed to initialize the random number generator used in MC eval
call MC_integration(n,3.0)

We need to initialize the random number generator by calling one of F90's built-in subroutines "call random_seed()". This built in random number (in reality, a pseudorandom number) generator is tolerable for the purposes of the class…better ones have been designed by computer scientists.

We then call the subroutine MC_integration. We first define a few variables, as you've seen before. One twist is that we want to keep extra bits of precision when we compute $f$ and $f^2$. To do this we can define these as

real (kind=8) :: f, f2

which doubles the internal representation of these by the computer.

After defining the variables, the heart of this subroutine is:

integral = 0.0

f= 0.0d0
f2 = 0.0d0

do i=1,n
    call random_number(x)
     x=x*end_val  ! random_number only returns uniformly distributed from [0.0, 1.0]

     f = f+integrand(x)
     f2 = f2+ (integrand(x)**2)
end do

where we first initialize our variables to zero (just to make sure!), and then perform a loop where:

1. Get a random number using the F90 built in subroutine "call random_number(x)". This puts a uniformly distributed number from $[0, 1]$ in the variable x.
2. We then have to multiply this to scale it over the base width of our integral
3. We then keep running sums of $f$ and $f^2$.

Finally, after repeating this for $n$ times, we evaluate our averages, and write out the value for the integral and its error:

f=f/n
f2=f2/n

integral = (end_val-0.0)*f
integral_err =  (end_val-0.0)*SQRT((f2 - f**2)/n)
write (*,*) "# MC integration = ",integral,"+/-",integral_err