316 Chapter 7.Random Numbers 7.8 Adaptive and Recursive Monte Carlo Methods This section discusses more advanced techniques of Monte Carlo integration.As examples of the use of these techniques,we include two rather different,fairly sophisticated, multidimensional Monte Carlo codes:vegas[1,2],and miser [3].The techniques that we discuss all fall under the general rubric of reduction of variance ($7.6),but are otherwise quite distinct. 三 Importance Sampling The use of importance sampling was already implicit in equations (7.6.6)and (7.6.7). We now return to it in a slightly more formal way.Suppose that an integrand f can be written as the product of a function h that is almost constant times another,positive,function g.Then its integral over a multidimensional volume V is fdv =(f/g)gdv hgdv (7.8.1) In equation(7.6.7)we interpreted equation (7.8.1)as suggesting a change of variable to G,the indefinite integral of g.That made gdV a perfect differential.We then proceeded to use the basic theorem of Monte Carlo integration,equation (7.6.1).A more general ad出 Press. THE interpretation of equation (7.8.1)is that we can integrate f by instead sampling h not, however,with uniform probability density dV,but rather with nonuniform density gdV.In this second interpretation,the first interpretation follows as the special case,where the means of generating the nonuniform sampling of gdv is via the transformation method,using the Programs indefinite integral G(see $7.2). More directly,one can go back and generalize the basic theorem (7.6.1)to the case of nonuniform sampling:Suppose that points z are chosen within the volume V with a probability density p satisfying 6 pdV=1 (7.8.2) 1C% The generalized fundamental theorem is that the integral of any function f is estimated,using N sample points 1,...,N,by f2/p2)-U/p)2 Numerical Recipes 10-621 I=fav= (7.8.3) uction. 43108 where angle brackets denote arithmetic means over the N points,exactly as in equation (7.6.2).As in equation (7.6.1),the "plus-or-minus"term is a one standard deviation error estimate.Notice that equation (7.6.1)is in fact the special case of equation (7.8.3),with (outside p constant 1/V. What is the best choice for the sampling density p?Intuitively,we have already seen Software. that the idea is to make h=f/p as close to constant as possible.We can be more rigorous by focusing on the numerator inside the square root in equation (7.8.3),which is the variance per sample point.Both angle brackets are themselves Monte Carlo estimators of integrals, so we can write s=(〉-(}≈∫sw-V5'=5w-ra (7.84) We now find the optimal p subject to the constraint equation (7.8.2)by the functiona ariation o-(Ew-ffrav]'+xfoa (7.8.5)316 Chapter 7. Random Numbers Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 7.8 Adaptive and Recursive Monte Carlo Methods This section discusses more advanced techniques of Monte Carlo integration. As examples of the use of these techniques, we include two rather different, fairly sophisticated, multidimensional Monte Carlo codes: vegas [1,2], and miser [3]. The techniques that we discuss all fall under the general rubric of reduction of variance (§7.6), but are otherwise quite distinct. Importance Sampling The use of importance sampling was already implicit in equations (7.6.6) and (7.6.7). We now return to it in a slightly more formal way. Suppose that an integrand f can be written as the product of a function h that is almost constant times another, positive, function g. Then its integral over a multidimensional volume V is
f dV =
(f /g) gdV =
h gdV (7.8.1) In equation (7.6.7) we interpreted equation (7.8.1) as suggesting a change of variable to G, the indefinite integral of g. That made gdV a perfect differential. We then proceeded to use the basic theorem of Monte Carlo integration, equation (7.6.1). A more general interpretation of equation (7.8.1) is that we can integrate f by instead sampling h — not, however, with uniform probability density dV , but rather with nonuniform density gdV . In this second interpretation, the first interpretation follows as the special case, where the means of generating the nonuniform sampling of gdV is via the transformation method, using the indefinite integral G (see §7.2). More directly, one can go back and generalize the basic theorem (7.6.1) to the case of nonuniform sampling: Suppose that points xi are chosen within the volume V with a probability density p satisfying
p dV =1 (7.8.2) The generalized fundamental theorem is that the integral of any function f is estimated, using N sample points x1,...,xN , by I ≡
f dV =
f p pdV ≈ f p ±  f 2/p2−f /p 2 N (7.8.3) where angle brackets denote arithmetic means over the N points, exactly as in equation (7.6.2). As in equation (7.6.1), the “plus-or-minus” term is a one standard deviation error estimate. Notice that equation (7.6.1) is in fact the special case of equation (7.8.3), with p = constant = 1/V . What is the best choice for the sampling density p? Intuitively, we have already seen that the idea is to make h = f /p as close to constant as possible. We can be more rigorous by focusing on the numerator inside the square root in equation (7.8.3), which is the variance per sample point. Both angle brackets are themselves Monte Carlo estimators of integrals, so we can write S ≡ f 2 p2 − f p 2 ≈
f 2 p2 pdV − 
f p pdV 2 =
f 2 p dV − 
f dV 2 (7.8.4) We now find the optimal p subject to the constraint equation (7.8.2) by the functional variation 0 = δ δp 
f 2 p dV − 
f dV 2 + λ
p dV  (7.8.5)