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)`

Here is one approach:

```
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))
```

Ouput:

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