Estimating the error function using Monte Carlo integration

Question

How do I perform Monte Carlo integration on the error function

$erf(x) = \frac{2}{\sqrt(\pi)}\int_{0}^{x}e^{-t^2}dt$

I have to estimate this using python but I'm not sure how the integral works at all. Can anyone help?

Answer 1

Monte Carlo integration is a technique used to estimate integrals. The key is to write the integrand as the product of a probability density function and some other function so that the integral effectively becomes an expected value problem. In your case:

$$ \operatorname{erf}(x) = \frac{2 x}{\sqrt\pi}\int_{0}^{x}e^{-t^2}/x\,\mathrm dt=\frac{2 x}{\sqrt\pi}\mathsf E(\exp(-X^2)), $$

where $X\sim\operatorname{Uniform}(0,x)$. Hence, we may estimate the error function by sampling $n$ psudeo-random variates from the $\operatorname{Uniform}(0,x)$ density and then computing $$ \widehat{\operatorname{erf}(x)}=\frac{2x}{n\sqrt\pi}\sum_{k=1}^n\exp(-x_i^2) $$

In MATLAB we may implement this as

data=rand(1e6,1);
erfEst=@(x) 2*x/sqrt(pi)*mean(exp(-(x*data).^2));