Numerical methods. Note 11
[~Home~] Monte Carlo Integration. Plain Monte Carlo. Numerical Methods. Note~ « 11 »

Random sampling. Plain Monte Carlo. We sample N points randomly from the integration region and evaluate the average of the integrand. The integral is equal to the volume of the integration region times the average of the integrand

ab  f(x) dx  ≈  V \mean{f} = V 1/N~i f(xi),
where V is the integration volume: V=b-a. For a n-dimentional rectangular region {x[] : a[i] <= x[i] <= b[i], i=1:n} the volume is V~=~∏(i=1: n)(bi-ai). The error estimate is (according to the central limit theorem) close to V~σ/√(N) where σ is the variance σ2~=~\mean{f2}~-~\mean{f}2.

Here is a typical algorithm for a n-dimensional plain Monte Carlo integration:

function plain_mc(n, a[n], b[n], f, N) returns [I,err] 
	for i=1..N :
		for k=1..n : x[i] = a[i] + random()*(b[i]-a[i])
		sum += f(x[])
		sum2 += f(x[])^2
	average = sum/N
	variance = sqrt( sum2/N - average^2 )
	vol = product( b[i]-a[i], i=1..n)
	I = vol*average
	err = vol*variance/sqrt(N)
The variance can be reduced by a cleverer sampling.

Importance sampling. Metropolis algorithm. Distribute points not uniformly but with larger density in the regions which contribute most (that is where the function is largest). Suppose that the points are distributed with the density ρ, that is the number of points Δn in the interval Δx is equal Δn = ρ(x)Δx. Then the estimate of the integral I is

I = ∑ fi Δxi = ∑ fi 1/ρi = ∑ fi/ρi
The corresponding variance is now σ2=\mean{(f/ρ)2}-\mean{f/ρ}2. Apparently if the ratio f/ρ is close to a constant the variance is reduced. It would be best to take ρ=|f| and sample directly from the function. However in practice it is typically expensive to evaluate the integrand. Therefore one builds the approximate density in a factorized form ρ(x,y,...,z) = ρx(x)ρy(y)...ρz(z). The standard VEGAS routine implements this variant of importance sampling.

Problems

  1. Implement a plain Monte Carlo multi-dimensional integration. The routine should calculate the average and the variance for a given function in multi dimensions over a given rectangular volume.
  2. Calculate 0π dx/π~0π dy/π~0π dz/π (1-cos(x)cos(y)cos(z))-1~=~Γ(1/4)4/(4π3) ≈  1.3932039296856768591842462603255

"Copyleft" © 2004 D.V.Fedorov (fedorov at phys dot au dot dk)