Homework "Monte Carlo integration"

1. (6 points) Plain Monte Carlo integration

   • Implement plain Monte Carlo multi-dimensional integration, e.g. along the following lines:
     static (double,double) plainmc(Func<vector,double> f,vector a,vector b,int N){ int dim=a.size; double V=1; for(int i=0;i<dim;i++)V*=b[i]-a[i]; double sum=0,sum2=0; var x=new vector(dim); var rnd=new Random(); for(int i=0;i<N;i++){ for(int k=0;k<dim;k++)x[k]=a[k]+rnd.NextDouble()*(b[k]-a[k]); double fx=f(x); sum+=fx; sum2+=fx*fx; } double mean=sum/N, sigma=Sqrt(sum2/N-mean*mean); var result=(mean*V,sigma*V/Sqrt(N)); return result; }
   • Calculate some interesting two-dimensional integrals with your Monte-Carlo routine (for example, calculate the area of a unit circle). Plot the estimated error and the actual error as functions of the number of sampling points. Check whether the actual error scales as ¹/_√N.
   • Try calculate

     ∫0π  dx/π ∫0π  dy/π ∫0π  dz/π [1-cos(x)cos(y)cos(z)]-1 = Γ(1/4)4/(4π3) ≈ 
     1.3932039296856768591842462603255
     

     This is a difficult singular integral, do not use it for your tests and error estimates, use something less sungular.

2. (3 points) Quasi-random sequences

   • Implement a multidimensional Monte-Carlo integrator that uses low-discrepancy (quasi-random) sequences. The error could be estimated by using two different sequences. Compare the scaling of the error with your pseudo-random Monte-Carlo integrator.

3. (1 point) Stratified sampling

   • Implement recursive stratified sampling, the variant where the subroutine receives the number of points N as an argument. Something like

     If N<nmin return N-point plain Monte Carlo estimate of integral and variance;
     Sample nmin points and estimate the integral and the variance;
     Find the dimension with largest sub-variance;
     Subdivide the volume along this dimension;
     Divide the remaning points between the two sub-volumes proportional to sub-variances;
     Dispatch two recursive calls on the sub-volumes;
     Estimate the grand integral and grand error;
     Return the grand integral and the grand error;