"Stochastic integration"
"Use NIntegrate to obtain the integral from 0 to 5 so that we can check the correctness of the code developed for stochasitc integration"
f[x_]:=x*(1+x^2)^(-3/2);
Integrate[f[x],x]
xlast=5.0;xinit=0.0;
pltf=Plot[f[x],{x,xinit,xlast},Filling->Axis];
Print[pltf]
nintegratef=NIntegrate[f[x],{x,xinit,xlast}];
Print[nintegratef]

(*First, build the rectangular box with dimension fmax x L *)
findmax=FindMaximum[f[x],{x,xlast}]
fmax=findmax[[1]];
L=xlast-xinit;
A=fmax*L;

(*In a for loop running from 1 to totalcountmax, generate a pair of random  real numbers (xrand,yrand) located inside the rectangular box. Add the 1 to inboxcount or otherwise. *)
totalcountmax=10000;

For[i=1,i≤totalcountmax,i++,
inbox[i]={0,0};
missed[i]={0,0};
];(*end for i*)


inboxcount=0;countmissed=0;
For[i=1,i≤totalcountmax,i++,
xrand=RandomReal[{xinit,xlast}];
yrand=RandomReal[{0,fmax}];
If[0<=yrand<=f[xrand],
inbox[i]={xrand,yrand};
inboxcount=inboxcount+1;
,(*else, record the missed points*)
missed[i]={xrand,yrand};
countmissed=countmissed+1;
];(*end if*)
Countinbox[i]=inboxcount;
Countmissed[i]=countmissed;
];(*end for i*)


tinbox[j_]:=Table[inbox[i],{i,1,j}];
tmissed[j_]:=Table[missed[i],{i,1,j}];
gp[j_]:=Graphics[{Blue,PointSize[0.01],Point[tinbox[j]]}];
gpm[j_]:=Graphics[{Red,PointSize[0.01],Point[tmissed[j]]}
];
pltfj[j_]:=Plot[f[x],{x,xinit,xlast},Frame->True,PlotLabel->{j,Style[Countmissed[j],Red],
Style[Countinbox[j],Blue]
}];

Manipulate[Show[pltfj[j],gpm[j],gp[j]],{j,1,totalcountmax,1}]
integratestochastic=A*inboxcount/totalcountmax;
Print["{nintegratef,integratestochastic,A,totalcountmax,inboxcount}",{nintegratef,integratestochastic,A,totalcountmax,inboxcount}];