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

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.

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