# Mathematica or math question: why my own Monte Carlo gives different result from the Mathematica's one?

Consider some function of arguments x,y,z, where the latter have some domain of definition $$\mathcal{R}$$:

function[x_, y_, z_] = Exp[-x^2]/(y^2 + 1.23)*Exp[-z];
zmin[x_, y_] = x + y*Exp[-y^2];
zmax[x_, y_] = 3.4*(x + 1.2*y*Exp[-y]);
ymin[x_] = 0.1*x*Exp[-x^3 - x^2] + 0.1*Exp[-Cos[x]^2];
ymax[x_] = x*Exp[-x^((2/3.)) + 1] + 0.22*Exp[-Cos[x]^4];
xmin = 0;
xmax = 10;


I need to evaluate an integral $$\mathcal{I} = \int \limits_{x_{\text{min}}}^{x_{\text{max}}}dx \ \int \limits_{y_{\text{min}}}^{y_{\text{max}}}dy \ \int \limits_{z_{\text{min}}}^{z_{\text{max}}}dz \ f(x,y,z)$$

For some reason, I would like to implement my own simple Monte Carlo integration. Namely, I start with the formula

$$\tag 1 \mathcal{I} \approx \frac{\Phi}{N}\sum_{i = 1}^{N}f(x_{i},y_{i},z_{i}), \quad \text{where} \quad \Phi = \int \limits_{x_{\text{min}}}^{x_{\text{max}}}dx\ \int \limits_{y_{\text{min}}}^{y_{\text{max}}}dy \ \int \limits_{z_{\text{min}}}^{z_{\text{max}}}dz, \quad (x_{i},y_{i},z_{i}) \in \mathcal{R}$$

First, I evaluate the phase space $$\Phi$$:

Measure =
NIntegrate[
1, {x, xmin, xmax}, {y, ymin[x], ymax[x]}, {z, zmin[x, y],
zmax[x, y]}]


84.0468

Next, I generate random points belonging to $$\mathcal{R}$$:

RandomxValues[Nevents_] := RandomReal[{xmin, xmax}, Nevents]
RandomyValues =
Hold@Compile[{{RandomxValuesTable, _Real, 1}},
Table[RandomReal[{ymin[RandomxValuesTable[[i]]],
ymax[RandomxValuesTable[[i]]]}], {i, 1,
Length[RandomxValuesTable], 1}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"] /. DownValues@ymin /.
DownValues@ymax // ReleaseHold;
RandomzValues =
Hold@Compile[{{RandomxValuesTable, _Real,
1}, {RandomyValuesTable, _Real, 1}},
Table[RandomReal[{zmin[RandomxValuesTable[[i]],
RandomyValuesTable[[i]]],
zmax[RandomxValuesTable[[i]], RandomyValuesTable[[i]]]}], {i,
1, Length[RandomxValuesTable], 1}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"] /. DownValues@zmin /.
DownValues@zmax // ReleaseHold;
RandomxyzPoints[Nevents_] := Block[{},
RandomxValuesTable = Quiet[RandomxValues[Nevents]];
RandomyValuesTable = Quiet[RandomyValues[RandomxValuesTable]];
RandomzValuesTable =
Quiet[RandomzValues[RandomxValuesTable, RandomyValuesTable]];
Join[Partition[RandomxValuesTable, 1],
Partition[RandomyValuesTable, 1], Partition[RandomzValuesTable, 1],
2]
]
(*An example demonstrating the performance of RandomxyzPoints*)
randomPointsTable = RandomxyzPoints[10^6]; // AbsoluteTiming


{0.209248,Null}

It looks like the generator works properly:

randomPointsTable1 = RandomxyzPoints[10^4];
Show[ListPointPlot3D[randomPointsTable1],
DiscretizeRegion[
ImplicitRegion[
xmin < x < xmax && ymin[x] < y < ymax[x] &&
zmin[x, y] < z < zmax[x, y], {x, y, z}](*,{{xmin,xmax},{0,30},{0,
30}}*)]]


Finally, my Monte Carlo integral is

FunctionSum =
Hold@Compile[{{randomPointsTable, _Real, 2}},
Sum[function[randomPointsTable[[i]][[1]],
randomPointsTable[[i]][[2]], randomPointsTable[[i]][[3]]], {i,
1, Length[randomPointsTable], 1}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"] /. DownValues@function // ReleaseHold
MonteCarloMy[Nevents_] :=
Measure/Nevents*FunctionSum[RandomxyzPoints[Nevents]]


Comparing with

MonteCarloMathematica =
NIntegrate[
function[x, y, z], {x, xmin, xmax}, {y, ymin[x], ymax[x]}, {z,
zmin[x, y], zmax[x, y]}, Method -> "MonteCarlo"]


0.131921

I find

MonteCarloMy[10^6]


1.49451

My question is: am I stupid since Eq. (1) is incorrect, or I made something wrong in the implementation?

An edit.

Let us simplify the domain of the definition of x,y,z, such that they are all independent:

zmin[x_, y_] = 0;
zmax[x_, y_] = 13;
ymin[x_] = 1.34;
ymax[x_] = 144;
xmin = 0;
xmax = 10;


and relaunch all the remaining code. The results, however, are still different:

MonteCarloMy[10^6]/MonteCarloMathematica

0.8007


I do not understand the reason for this discrepancy.

Edit 2.

So, if continuing with my method, following comments and the answer, I need to define $$\Delta z[x,y] = \int\limits_{\mathcal{R}_{x,y}} dz$$, $$\Delta y[x] = \int\limits_{\mathcal{R}_{x}} dy$$, and instead of pure sum of function I need to multiply it with the weights $$\Delta z \times \Delta y$$, and divide by the sum of $$\Delta z \times \Delta y$$.

• I think the simple answer is that your x, y, z values are not drawn independently. That is, y depends on x and z depends on both x and y. If you would draw them independently, then an x-value should be very rare if its corresponding y-z-slice is small. But in your code all x-values in the interval [xmin, xmax] seem to have the same probability. Aug 7, 2022 at 19:43
• @HenrikSchumacher Thanks! Do you mean that I need to compute the weight $w(x_{i}) = \int dydz$, and then re-generate the points distribution according to the value of weight? Aug 7, 2022 at 19:48
• Yes, that might work. In fact, you parameterize your domain by a cuboid. So what you basically have to do is figure out this mapping and reweigh your sampling the mapping's Jacobian determinant. Aug 7, 2022 at 19:52
• @HenrikSchumacher : how do you think, should it work if instead of Eq. (1) I would write $$\mathcal{I} = \frac{\Phi}{N} \times \frac{\sum_{i}\omega_{x_{i}}f(x_{i},y_{i},z_{i})}{\sum_{i}\omega_{x_{i}}},$$ where the weights $\omega_{x_{i}}$ are $$\omega_{x_{i}} = \int dx \ dy \quad \text{at} \quad x_{i} ?$$ Aug 7, 2022 at 20:06
• You could sample (x,y,z) uniformly from a cuboid that contains your region, then reject the sample if it lies outside of your region. This will give you the region measure and the function integral at the same time. en.wikipedia.org/wiki/Rejection_sampling Aug 7, 2022 at 20:12

There are two issues at play: Non-uniform sampling and bad convergence.

First, let's consider your example with independent region bounds:

zmin[x_, y_] = 0;
zmax[x_, y_] = 13;
ymin[x_] = 1.34;
ymax[x_] = 144;
xmin = 0;
xmax = 10;

MonteCarloMathematica =
NIntegrate[
function[x, y, z], {x, xmin, xmax}, {y, ymin[x], ymax[x]}, {z,
zmin[x, y], zmax[x, y]}, Method -> "MonteCarlo"]
(* NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 0.4398652196837818 and 0.06682921207520477 for the integral and error estimates. *)
(* 0.439865 *)

MonteCarloMy[10^5]
(* 0.509626 *)


As you note, the results don't appear to match. However, read the warning produced by Mathematica: It estimates that the integral (including error) is 0.440 +- 0.067, which actually matches your result. This is a good point to emphasize again that you should not ignore warning messages produced by Mathematica, especially not numerical methods that might want to warn you about issues. Giving both methods more sampling points reduces the discrepancy:

MonteCarloMathematica =
NIntegrate[
function[x, y, z], {x, xmin, xmax}, {y, ymin[x], ymax[x]}, {z,
zmin[x, y], zmax[x, y]}, Method -> "MonteCarlo", MaxPoints -> 10^8]
(* 0.548626 *)

MonteCarloMy[10^8]
(* 0.546474 *)


Next, to your original situation: As noted by others, your sampling is far from uniform:

function[x_, y_, z_] = Exp[-x^2]/(y^2 + 1.23)*Exp[-z];
zmin[x_, y_] = x + y*Exp[-y^2];
zmax[x_, y_] = 3.4*(x + 1.2*y*Exp[-y]);
ymin[x_] = 0.1*x*Exp[-x^3 - x^2] + 0.1*Exp[-Cos[x]^2];
ymax[x_] = x*Exp[-x^((2/3.)) + 1] + 0.22*Exp[-Cos[x]^4];
xmin = 0;
xmax = 10;

Graphics3D[{AbsolutePointSize@1, Point@RandomxyzPoints[10^5]},
BoxRatios -> {1, 1, 1}]


Using (highly unoptimized) rejection sampling, we get the following:

function[x_, y_, z_] = Exp[-x^2]/(y^2 + 1.23)*Exp[-z];
zmin[x_, y_] = x + y*Exp[-y^2];
zmax[x_, y_] = 3.4*(x + 1.2*y*Exp[-y]);
ymin[x_] = 0.1*x*Exp[-x^3 - x^2] + 0.1*Exp[-Cos[x]^2];
ymax[x_] = x*Exp[-x^((2/3.)) + 1] + 0.22*Exp[-Cos[x]^4];
xmin = 0;
xmax = 10;

yminGlobal = First@Minimize[{ymin[x], xmin <= x <= xmax}, x]
(* 0.0367879 *)

ymaxGlobal = First@Maximize[{ymax[x], xmin <= x <= xmax}, x]
(* 1.33353 *)

zminGlobal =
First@Minimize[{zmin[x, y], xmin <= x <= xmax,
yminGlobal <= y <= ymaxGlobal}, {x, y}]
(* 0.0367382 *)

zmaxGlobal =
First@Maximize[{zmax[x, y], xmin <= x <= xmax,
yminGlobal <= y <= ymaxGlobal}, {x, y}]
(* 35.5009 *)

RandomxyzPoint :=
With[{x = RandomReal[{xmin, xmax}],
y = RandomReal[{yminGlobal, ymaxGlobal}],
z = RandomReal[{zminGlobal, zmaxGlobal}]},
Quiet@If[
ymin[x] <= y <= ymax[x] && zmin[x, y] <= z <= zmax[x, y], {x, y,
z}, RandomxyzPoint]
]

RandomxyzPoints[Nevents_] := Table[RandomxyzPoint, Nevents]

MonteCarloMathematica =
NIntegrate[
function[x, y, z], {x, xmin, xmax}, {y, ymin[x], ymax[x]}, {z,
zmin[x, y], zmax[x, y]}, Method -> "MonteCarlo", MaxPoints -> 10^6]
(* 0.132109 *)

MonteCarloMy[10^6]
(* 0.13359 *)

• Thanks! But in the case of non-uniform distribution, how did you get 0.13359 for MonteCarloMy? Without modifying the code, I get 1.49451. Aug 7, 2022 at 21:01
• @JohnTaylor Sorry, I forgot to actually paste the rejection sampling code into my answer... Should be fixed now Aug 7, 2022 at 21:12
• Thanks! This is an alternative to my approach modified by the second edit (Edit 2). Aug 7, 2022 at 21:19