Validating results of a Monte Carlo integration

Question

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

Answer 1

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

Answer 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.

Answer 3

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.