Wolfram Mathematica Monte Carlo for integrals approximation

Question

I wanted to implement the Monte Carlo method for multiple integrals approximation in Wolfram Mathematica. Namely I wanted to let the user insert as input the dimension of the integral and the number of samples to be considerated and then perform the method. I implemented the following code:

monteCarloIntegration[f_, a_, b_, n_, dim_] :=
Module[{randomNumbers, values, average, integrand},
randomNumbers = RandomReal[#, n] & /@ Transpose[{a, b}^dim];
integrand = Apply[f, randomNumbers, {1}];
values = Total[integrand];
average = Mean[values];
(b - a)^dim*average];

(Example for dim = 3 (triple integral)):

f[x_, y_, z_] := x^2 + y^2 + z^2;
a = {0, 0, 0};
b = {1, 1, 1};
n = 1000;
dim = 3;
trueV = (1/4) (Pi)^(3/2);
(*monteCarloIntegration[f,a,b,n,dim]*)

But when I plot the following:

 errorValues = Table[Abs[trueV-monteCarloIntegration[f, a, b, i, dim]], {i, 1, 1000}];
 ListLinePlot[errorValues,
 PlotLabel -> "Error between the approximation and the real value",
 AxesLabel -> {"Number of Samples", "Error"}]

the graphic is displayed (I mean the axis and the label) but it's an empty plot. Could someone hlep me to fix this? Thank u.

Answer 1

The problem is that your f is being applied to the entire list of random values, which varies in size based on "n", which is set to "i" in the Table. So you apply f to a single random value for i=1, and to lists of length 1000 for i=1000. Since f only accepts 3 arguments, it remains unevaluated most of the time. How do you want to handle this? I.e. generate 1000 sets of 3, partition each list into sets of three and deal with leftovers? Etc. The built-in NIntegrate does include a Monte Carlo method by the way, which should handle your integrand.

I'm going to assume that you want to generate i=1,1000 groups of 3. If I replace your monteCarloIntegration with this:

monteCarloIntegration[f_, a_, b_, n_, dim_] := 
  Module[{randomNumbers, values, average, integrand}, 
   randomNumbers = RandomReal[#, {n, 3}] & /@ Transpose[{a, b}^dim];
   integrand = Apply[f, randomNumbers, {2}];
   values = Total[integrand];
   average = Mean[values];
   (b - a)^dim*average];

I get the following result for the plot:

However, I believe there's another problem with your approach. The basic idea of the Monte Carlo method is to enclose the region of interest inside a larger region which has an easy to calculate area. You then check to see if a given random point is inside the boundary of the difficult to compute region, and take the ratio of the number of points that land inside to the total number of random trials. This equals the ratio of areas inside/outside. You don't seem to be doing this at present.

As I mentioned, Monte Carlo methods (MonteCarlo, AdaptiveMonteCarlo,QuasiMonteCarlo, AdaptiveQuasiMonteCarlo) are supported by NIntegrate. For example,

NIntegrate[
   Boole[x^2 + y^2 + z^2 < 1], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, 
   Method -> #] & /@ {"MonteCarlo", "AdaptiveMonteCarlo", 
  "QuasiMonteCarlo", "AdaptiveQuasiMonteCarlo"}

Details are in the documentation.

Incidentally,your "trueV" seems to be off. I get Pi/6 for the integral instead of 1/4 (Pi)^(3/2):

Integrate[Boole[x^2 + y^2 + z^2 < 1], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]