Numerical methods. Note 13


	[ Home ]	Numerical Methods. Note « 13 »

NB: There have been suggested two subjects for the last lecture:

Two-point boundary value problem.
Monte-Carlo methods in physics.

Monte Carlo integration.

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

∫_{_a}^{^b} f(x) dx ≈ (b-a)<f> = (b-a) ¹/_N ∑_{_i} f(x_i) .

The error estimate δ is (probably) close to (b-a)^σ/_√(N) where σ is the variance σ²=<f²>-<f>². 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 = ∑ f_i Δx_i = ∑ f_i ¹/_{ρ_i} = ∑ ^f_i/_{ρ_i}

The corresponding variance is now σ²=<(f/ρ)²>-<f/ρ>². 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.

Stratified sampling. Distribute more points in the regions where the variance is largest. The MISER routine implements recursive stratified sampling along the following lines

sample few points randomly from the given volume
calculate average and variance
if variance is small enough then
	estimate the grand average and grand variance; done
else
	estimate which dimension gives the largest variance
	subdivide your volume in two along this dimension
	recurse for each of the sub-volumes

Problems

Obligatory exercise: Implement a plain Monte Carlo multi-dimensional integration. The routine should calculate an integral of a function in multi dimensions over a given rectangular volume.
Non-obligatory exercise: stratified sampling.