Monte Carlo Integration

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Straightforward implementation without the weight function

Simple Monte Carlo method approximates a multivariable integral by

I=(∫_Ω)f((Subscript[x, 1]),…,(Subscript[x, d]))(d^d)x=((|Ω|)/N)(∑_(i=1)^N)f((Subscript[x, 1]^((i))),…,(Subscript[x, d]^((i))))

where | Ω | denotes the volume of the region Ω.
The following function approximates a dim-dimensional integral of a function f over the hypercube with boundaries given by dim-dimensional vectors a and b. It returns as an output an array with 3 rows (the number of points, approximation of the integral, error estimation). The number of columns is exp, while the number of points for each approximation is given by  N=(base^exp). The error estimation is obtain as

probable error=0.67((|Ω|)/(Sqrt[N]))((((Underoverscript[(1/N)∑, (i=1), N])f((((Subscript[x, 1]^((i))),…,(Subscript[x, d]^((i)))))^2)-((((Underoverscript[(1/N)∑, (i=1), N])f((Subscript[x, 1]^((i))),…,(Subscript[x, d]^((i))))))^2)))^(1/2))

mySimpleMCIntND[dim_,f_,a_,b_,base_,exp_]:=
Module[{i,ie,rndPt,volume,I,II,ip,n,out},
out=ConstantArray[0,{3,exp}];
volume=Times@@(b-a);
n=base^exp; (* the largest number of points *)
I=0.0;II=0.0;
ie=1; (* exponent, only results for base^ie points are returned *)
Do[
rndPt=a+(b-a)RandomReal[{0.0,1.0},dim];
I=I+f[rndPt]; (* sum of function values *)
II=II+f[rndPt]^2;(* sum of squares of function values *)
If[Mod[i,base^ie]==0,
out[[All,ie]]={i,volume*I/i,0.67*volume*Sqrt[(II/i-(((I/i))^2))/i]};
ie++;
],
{i,1,n}
];
Return[out];
];

The volume of the unit ball in dim dimensions:

(* the function to integrate, Norm[x] returs the length of vector x *)
f1[x_]:=If[Norm[x]≤1.0,1.0,0.0]; 
dim=2; (* dimension *)
a=ConstantArray[-1,{dim}];b=-a;(* boundaries, hypercube arround the origin *)
int=mySimpleMCIntND[dim,f1,a,b,4,9];(* the largest number of points is 4^10 *)
plot=ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-(π^(dim/2))/Gamma[1+dim/2]]}]},
AxesLabel->{"N","absolute error"},PlotLegends->{"Estimated Error","Actual Error"}];
Print["#points integral   error"];
Print[MatrixForm[Transpose[int]],"   ",Show[plot,ImageSize->{300,180}]];
Print["Exact values: hypercube volume = ",N[2^dim],", unit ball volume = ",N[(π^(dim/2))/Gamma[1+dim/2]]];

Integration of the function 1/(1+x^2)

f[x_]:=1/(1+Norm[x]^2);  
dim=1;
a=ConstantArray[0,{dim}];b=ConstantArray[1,{dim}];
int=mySimpleMCIntND[dim,f,a,b,4,10];
plot=ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-π/4]}]},AxesLabel->{"N","absolute error"},PlotLegends->{"Estimated Error","Actual Error"}];
Print["#points integral   error"];
Print[MatrixForm[Transpose[int]],"   ",Show[plot,ImageSize->{300,180}]];
Print["Exact: ",N[π/4,10],", Error: ",N[Abs[int[[2,10]]-π/4],10]];

Integration of the function 1/(1+x^2) with the linear function as a weight function

f[x_]:=1/(1+Norm[x]^2);
w[x_]:=(4-2Norm[x])/3; (* normalized linear function used as a weight function *)
fx[y_]:=2-Sqrt[4-3Norm[y]]; (* inverse of the CDF of the density w(x) *)
f1[y_]:=f[fx[y]]/w[fx[y]]; (* actual function to integrate *)
Plot[{f[{x}],w[{x}],f[{x}]/w[{x}]},{x,a[[1]],b[[1]]},PlotRange->{0,1.5},AxesLabel->{"x","f(x)"},PlotLegends->{"f(x)","w(x)","f(x)/w(x)"}]
a={0};b={1};
dim=1;
int = mySimpleMCIntND[dim,f1,a,b,4,10];
errorPlot=ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-π/4]}]},AxesLabel->{"N","absolute error"},PlotLegends->{"Estimated Error","Actual Error"}];
Print["#points integral   error"];
Print[MatrixForm[Transpose[int]],"   ",Show[errorPlot,ImageSize->{300,180}]];
Print["Exact: ",N[π/4,10],", Error: ",N[Abs[int[[2,10]]-π/4],10]];