# Validating results of a Monte Carlo integration

Suppose I want to use the Monte Carlo integration method to compute the following integral

$\int_D (e^{x^{2}} + e^{y^{2}}) \; dx \; dy$ where $D$ is some regular hexagon. I have managed to write the code and everything and the results I'm getting seem valid.

The problem is now that I do not know how can I show that the algorithm does what it's supposed to do. Is there another way to compute that integral? Is there any way to validate my results?

EDIT: I'm talking about the "sample-mean" method

-
This is just a guess, but maybe you have a typo and you really want to integrate over $e^{-(x^2 + y^2)}$ instead? This would be a 2d gaussian which is probably more common than the function you wrote. – opt Jan 25 '12 at 18:54
Perhaps see http://math.stackexchange.com/questions/1635250/ for more on Monte Carlo integration. – BruceET Feb 1 at 17:42

## 3 Answers

Let $f(x,y) = \mathrm{e}^{x^2} + \mathrm{e}^{y^2}$, and suppose we would like to integrate this function over a valid triangle with vertices at $p_1 = (x_1,y_1)$, $p_2 (x_2,y_2)$ and $p_3 = (x_3,y_3)$.

This can be computed in closed form, because $\int_D f(x,y) \mathrm{d} x \mathrm{d} y = \int_D \mathrm{e}^{x^2} \mathrm{d} y \mathrm{d} x + \int_D \mathrm{e}^{y^2} \mathrm{d} x \mathrm{d} y$ Now integration over $y$ for each fixed $x$ in the first integral, and over $x$ for each fixed $y$ in the second can be easily done, and will yield a linear function of $x$ and $y$ respectively.

Then you are down to $$\int_a^b (c x+d) \mathrm{e}^{x^2} \mathrm{d} x = \frac{c}{2} \left( \mathrm{e}^{b^2}- \mathrm{e}^{a^2} \right) + \frac{d \sqrt{\pi}}{2} \left( \operatorname{erfi}(b) - \operatorname{erfi}(a) \right)$$

Coming back to the setting stated at the beginning of the post, here is a little code in Mathematica that will do the computation exactly:

In[6]:= pred[{x_, y_}, {x1_, y1_}, {x2_, y2_}, {x3_,
y3_}] := ((x - x1) (y2 - y1) - (y - y1) (x2 - x1)) ((x3 - x1) (y2 -
y1) - (y3 - y1) (x2 - x1)) > 0

In[7]:= InTriangle[px_, p1_, p2_, p3_] :=
pred[px, p1, p2, p3] && pred[px, p2, p3, p1] && pred[px, p3, p1, p2]

In[10]:= int[p1 : {x1_, y1_}, p2 : {x2_, y2_}, p3 : {x3_, y3_}] :=
Block[{x, y},
Integrate[(Exp[x^2] + Exp[y^2]) Boole[
InTriangle[{x, y}, p1, p2, p3]], {x, Min[x1, x2, x3],
Max[x1, x2, x3]}, {y, Min[y1, y2, y3], Max[y1, y2, y3]}]]

In[13]:= int[{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}] // FunctionExpand

Out[13]= 1/6 (-Sqrt[3] + 6 Sqrt[3] E^(1/4) - 2 Sqrt[3] E^(3/4) -
3 Sqrt[3] E - 3 Sqrt[3 Pi ] Erfi[1/2] +
3 Sqrt[3 Pi] Erfi[1] + 3 Sqrt[ Pi] Erfi[Sqrt[3]/2])

-

The simplest way to validate your Monte Carlo results would be to use numerical quadrature of modestly accurate order on the region D. Since D is a regular hexagon (but otherwise unspecified), it can be split into six equilateral triangles.

Quadrature rules in one-dimension are well introduced in undergraduate numerical methods courses, but if you want, take a formula from this article on quadrature rules for triangles.

Typically a series of such approximations is done, either with increasingly small subdivisions of area or with increasingly higher-order accuracy formulas, so that the convergence of the approximations can be seen.

-

Because the integrand is positive you can easily find lower and upper bounds by integrating over nicer domains that are subsets or supersets of the hexagon domain.

-
While this is true, such an approach require fine mesh triangulation of the hexagon to achieve modest precision, as $\exp(x^2) + \exp(y^2)$ grow rather fast. – Sasha Jan 25 '12 at 18:40