10
$\begingroup$

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

enter image description here

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

Improve this question
$\endgroup$
7
  • 6
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Commented Aug 7, 2022 at 19:43
  • $\begingroup$ @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? $\endgroup$
    – John Taylor
    Commented Aug 7, 2022 at 19:48
  • 1
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Commented Aug 7, 2022 at 19:52
  • $\begingroup$ @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} ? $$ $\endgroup$
    – John Taylor
    Commented Aug 7, 2022 at 20:06
  • 5
    $\begingroup$ 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 $\endgroup$
    – Roman
    Commented Aug 7, 2022 at 20:12

1 Answer 1

Reset to default
10
$\begingroup$

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

enter image description here

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 *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – John Taylor
    Commented Aug 7, 2022 at 21:01
  • 1
    $\begingroup$ @JohnTaylor Sorry, I forgot to actually paste the rejection sampling code into my answer... Should be fixed now $\endgroup$
    – Lukas Lang
    Commented Aug 7, 2022 at 21:12
  • $\begingroup$ Thanks! This is an alternative to my approach modified by the second edit (Edit 2). $\endgroup$
    – John Taylor
    Commented Aug 7, 2022 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.