Monte Carlo Integration

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Straightforward implementation without the weight function

mySimpleMCIntND[dim_,f_,a_,b_,base_,exp_]:=
Module[{i,ie,rndPt,volume,I,II,ip,n=base^exp,out},
out=ConstantArray[0,{3,exp}];
volume=Times@@(b-a);
I=0;II=0;
ie=1;
Do[
rndPt=a+(b-a)RandomReal[{0.0,1.0},dim];
I=I+f[rndPt];
II=II+((f[rndPt])^2);
If[Mod[i,base^ie]==0,
out[[All,ie]]={i,volume I/i,volume Sqrt[(II/i-(((I/i))^2))/i]};
ie++;
],
{i,1,n}
];
Return[out];
];

Unit ball:

f1[x_]:=If[Norm[x]≤1.0,1.0,0.0]; 
dim=10;
a=ConstantArray[-1,{dim}];b=-a;
int=mySimpleMCIntND[dim,f1,a,b,4,10];
MatrixForm[Transpose[int]]
Print["Exact: ",N[{2^dim,(π^(dim/2))/Gamma[1+dim/2]},10]];

ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-(π^(dim/2))/Gamma[1+dim/2]]}]}]

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];
MatrixForm[Transpose[int]]
Print["Exact: ",N[π/4,10]];

ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-π/4]}]}]

Analytic weight function

f[x_]:=1/(1+Norm[x]^2);
w[x_]:=(4-2Norm[x])/3;
fx[y_]:=2-Sqrt[4-3Norm[y]];
f1[y_]:=f[fx[y]]/w[fx[y]];
a={0};b={1};
dim=1;
int = mySimpleMCIntND[dim,f1,a,b,4,10];
MatrixForm[Transpose[int]]
Print["Exact: ",N[π/4,10]];

ListLogLogPlot[{Transpose[{int[[1,All]],int[[3,All]]}],Transpose[{int[[1,All]],Abs[int[[2,All]]-π/4]}]}]

Metropolis-Hastings algorithm

Metropolis-Hastings algorithm for generating random numbers with a given probability density function ρ(x) starting in (Subscript[x, 0]) with a random step in a 2δ-cube around the current position:

myRandomMetropolisND[n_,dim_,ρ_,a_,b_,x0_,δ_]:=
Module[{x=x0,xtrial,new,ratio,ρx,out=ConstantArray[0.0,{n,dim}]},
ρx=ρ[x];
Do[
xtrial=x+δ RandomReal[{-1,1},dim];
new=1;
Do[
If[xtrial[[d]]<a[[d]]||xtrial[[d]]>b[[d]],new=0],
{d,1,dim}
];
If[new==1,
ratio=ρ[xtrial]/ρx;
If[ratio⩾1||ratio>RandomReal[],x=xtrial;ρx=ρx ratio]
];
out[[i]]=x,
{i,1,n}
];
Return[out]
];

Implementation with the weight function using Metropoli-Hastings algorithm

Monte Carlo integration of the function f((OverVector[x])) with a weight function ρ((OverVector[x]))

myMCInt[n_,dim_,f_,a_,b_,ρ_]:=Module[{x0,δ,rndPts,I,II},
x0=(a+b)/2;
δ=(b-a)/10;
rndPts=myRandomMetropolisND[n,dim,ρ,a,b,x0,δ];
I=(Plus@@(f[#]/ρ[#]&/@rndPts))/n;
II=(Plus@@((((f[#]/ρ[#]))^2)&/@rndPts))/n;
Return[{I,Sqrt[(II-(I^2))/n]}]
];

1D example

dim=1;
f[x_]:=1/(1+Norm[x]^2);
ρ[x_]:=(4-2Norm[x])/3;
n=10000;
a=ConstantArray[0,{dim}];b=ConstantArray[1,{dim}];
Int=N[myMCInt[n,dim,f,a,b,ρ],10]
Print["Exact: ",N[π/4,10], ", diff = ",N[Int[[1]]-π/4]];