Monte-Carlo Integration: the Error Function

Question

The error function, $\mathrm{erf}(x)$ is defined as $$ \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x\mathrm{e}^{-t^2}\;\mathrm{d}t. $$ The integral cannot be evaluated in closed form, but many numerical approximations and converging series for $\mathrm{erf}(x)$ are known. Python's math module provides the method erf() for its evaluation.

One approximate method for evaluating the error function is Monte Carlo integration. $f(t) = \mathrm{e}^{-t^2}$ is a monotonically decreasing function of $x$, and $\mathrm{e}^{0} = 1$ so consider the rectangular region of the $tf$-plane bounded by $0 \le t \le x$ and $f(x) \le f \le 1$. The area under curve $f(t)$ can be approximated by picking a large number of random points inside this rectangle, and counting how many of them lie under the curve as a proportion of the total number picked.

Write a program to estimate the value of $\mathrm{erf}(\frac{1}{2})$ and compare it with the result of the builtin math.erf(0.5)

Solution

import math
import random

N = 10000

def f(t):
    return math.exp(-t**2)

x = 0.5
print('math built-in erf(1/2) = {}'.format(math.erf(x)))

i = 0
I = 0.
ymin = f(x)
while i < N:
    i += 1
    xi, yi = x * random.random(), ymin + (1. - ymin) * random.random()
    if yi <= f(xi):
        I += 1.
I /= N
I *= x * (1. - ymin)
I += x * ymin
I *= 2./math.sqrt(math.pi)
print('Monte Carlo (N = {}) estimate for erf(1/2) = {}'.format(N, I))

math built-in erf(1/2) = 0.5204998778130465
Monte Carlo (N = 10000) estimate for erf(1/2) = 0.5197364311964482