f7-6

7.6SimpleMonteCarloIntegration 295Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Validating the Correctness of Hardware Implementations of the NBS Data Encryption Stan- dard,1980,NBSSpecialPublication500–20(Washington:U.S.DepartmentofCommerce, National Bureau of Standards). [3] Meyer,C.H.andMatyas,S.M.1982, Cryptography:ANewDimensioninComputerDataSecurity (New York: Wiley). [4] Knuth,D.E. 1973, SortingandSearching ,v ol.3of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), Chapter 6. [5] Vitter,J.S.,andChen,W-C.1987, DesignandAnalysisofCoalescedHashing (NewYork:Oxford University Press). [6] 7.6 Simple Monte Carlo Integration Inspirationsfornumericalmethodscanspringfromunlikelysources. “Splines” first were flexible strips of wood used by draftsmen. “Simulated annealing” (we shall see in §10.9)is rooted in a thermodynamicanalogy. And who does not feel at least a faint echo of glamor in the name “Monte Carlo method”? Supposethatwe pick Nrandompoints,uniformlydistributedina multidimen- sional volume V. Call them x1,...,x N. Then the basic theorem of Monte Carlo integrationestimates the integralofafunction foverthe multidimensionalvolume, /integraldisplay fd V≈V/angbracketleftf/angbracketright±V/radicalBigg /angbracketleftf2/angbracketright−/angbracketleftf/angbracketright2 N(7.6.1 ) Heretheanglebracketsdenotetakingthearithmeticmeanoverthe Nsamplepoints, /angbracketleftf/angbracketright≡1 NN/summationdisplay i=1f(xi)/angbracketleftbig f2/angbracketrightbig ≡1 NN/summationdisplay i=1f2(xi)( 7.6.2 ) The “plus-or-minus” term in (7.6.1) is a one standard deviation error estimate for the integral, not a rigorous bound; further, there is no guarantee that the error is distributedasaGaussian,sotheerrortermshouldbetakenonlyasaroughindication of probable error. Suppose that you want to integrate a function gover a region Wthat is not easy to sample randomly. For example, Wmight have a very complicated shape. No problem. Just find a region Vthatincludes Wand thatcaneasily be sampled (Figure 7.6.1), and then define fto be equal to gfor points in Wand equal to zero for points outside of W(but still inside the sampled V). You want to try to make Venclose Was closely as possible, because the zero values of fwill increase the error estimate term of (7.6.1). And well they should: points chosen outside of W have no information content, so the effective value of N, the number of points, is reduced. The error estimate in (7.6.1) takes this into account. General purpose routines for Monte Carlo integration are quite complicated (see§7.8),butaworkedexamplewill showtheunderlyingsimplicityofthemethod. Supposethatwe want to findthe weightandthe positionof thecenterof mass ofan 296 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).area A ∫fdx Figure 7.6.1. Monte Carlo integration. Random points are chosen within the area A. The integral of the function fis estimated as the area of Amultiplied by the fraction of random points that fall below the curve f.R efinements on this procedure can improve the accuracy of the method; see text. 02 424y x 1 Figure 7.6.2. Example of Monte Carlo integration (see text). The region of interest is a piece of a torus, bounded bythe intersection oftwo planes. Thelimits ofintegration ofthe region cannot easily bewritten in analytically closed form, so Monte Carlo is a useful technique. 7.6SimpleMonteCarloIntegration 297Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).object of complicated shape, namely the intersection of a torus with the edge of a largebox. Inparticularlettheobjectbede finedbythethreesimultaneousconditions z2+/parenleftBig/radicalbig x2+y2−3/parenrightBig2 ≤1( 7.6.3 ) (torus centered on the origin with major radius =4, minor radius =2) x≥1 y≥− 3( 7.6.4 ) (two faces of the box, see Figure 7.6.2). Suppose for the moment that the object has a constant density ρ. We wanttoestimatethefollowingintegralsovertheinteriorofthecomplicated object: /integraldisplay ρd xd yd z/integraldisplay xρ dx dy dz/integraldisplay yρ dx dy dz/integraldisplay zρ dx dy dz (7.6.5 ) The coordinates of the center of mass will be the ratio of the latter three integrals (linear moments) to the first one (the weight). In the followingfragment,the region V, enclosingthe piece-of-torus W, is the rectangular box extending from 1 to 4 in x,−3t o4i n y, and−1t o1i n z. n= Set to the number of sample points desired. den= Set to the constant value of the density. sw=0. Zero the various sums to be accumulated. swx=0.swy=0. swz=0. varw=0.varx=0.vary=0. varz=0. vol=3.*7.*2. Volume of the sampled region. do 11j=1,n x=1.+3.*ran2(idum) Pick a point randomly in the sampled region. y=-3.+7.*ran2(idum) z=-1.+2.*ran2(idum)if (z**2+(sqrt(x**2+y**2)-3.)**2.le.1.)then Is it in the torus? sw=sw+den If so, add to the various cumulants. swx=swx+x*denswy=swy+y*denswz=swz+z*den varw=varw+den**2 varx=varx+(x*den)**2vary=vary+(y*den)**2 varz=varz+(z*den)**2 endif enddo 11 w=vol*sw/n The values of the integrals (7.6.5), x=vol*swx/n y=vol*swy/nz=vol*swz/n dw=vol*sqrt((varw/n-(sw/n)**2)/n) and their corresponding error estimates. dx=vol*sqrt((varx/n-(swx/n)**2)/n)dy=vol*sqrt((vary/n-(swy/n)**2)/n)dz=vol*sqrt((varz/n-(swz/n)**2)/n) 298 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).A change of variable can often be extremely worthwhile in Monte Carlo integration. Suppose, for example, that we want to evaluate the same integrals,but for a piece-of-torus whose density is a strong function of z, in fact varying according to ρ(x, y, z )=e 5z(7.6.6 ) One way to do this is to put the statement den=exp(5.*z) inside the if...then block, just before denisfirst used. This will work, but it is a poor way to proceed. Since (7.6.6) falls so rapidly to zero as zdecreases (down to its lower limit −1), most sampled points contribute almost nothing to the sum of the weight or moments. These points are effectively wasted, almost as badly as those that fall outside of the region W. A change of variable, exactly as in the transformation methods of §7.2, solves this problem. Let ds=e5zdzsothat s=1 5e5z,z =1 5ln(5s)( 7.6.7 ) Then ρdz =ds, and the limits −1<z< 1become .00135 <s< 29.682. The program fragment now looks like this n= Set to the number of sample points desired. sw=0.swx=0. swy=0. swz=0.varw=0. varx=0. vary=0.varz=0.ss=(0.2*(exp(5.)-exp(-5.))) Interval of sto be random sampled. vol=3.*7.*ss Volume in x,y,s -space. do 11j=1,n x=1.+3.*ran2(idum)y=-3.+7.*ran2(idum) s=.00135+ss*ran2(idum) Pick a point in s. z=0.2*log(5.*s) Equation (7.6.7). if (z**2+(sqrt(x**2+y**2)-3.)**2.lt.1.)then sw=sw+1. Density is 1, since absorbed into definition of s. swx=swx+xswy=swy+yswz=swz+z varw=varw+1. varx=varx+x**2vary=vary+y**2 varz=varz+z**2 endif enddo 11 w=vol*sw/n The values of the integrals (7.6.5), x=vol*swx/n y=vol*swy/nz=vol*swz/n dw=vol*sqrt((varw/n-(sw/n)**2)/n) and their corresponding error estimates. dx=vol*sqrt((varx/n-(swx/n)**2)/n)dy=vol*sqrt((vary/n-(swy/n)**2)/n)dz=vol*sqrt((varz/n-(swz/n)**2)/n) 7.7Quasi-(thatis,Sub-)RandomSequences 299Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If youthinkfor a minute,youwill realize that equation(7.6.7)was useful only becausethepartoftheintegrandthatwewantedtoeliminate( e5z)wasbothintegrable analytically,andhadanintegralthatcouldbeanalyticallyinverted. (Compare §7.2.) In general these properties will not hold. Question: What then? Answer: Pull out of the integrand the “best”factor that canbe integrated and inverted. The criterion for“best”is to try to reduce the remaining integrand to a function that is as close as possible to constant. Thelimitingcase is instructive: Ifyoumanageto makethe integrand fexactly constant, and if the region V, of known volume, exactlyencloses the desired region W,thentheaverageof fthatyoucomputewillbeexactlyitsconstantvalue,andthe error estimate in equation (7.6.1) will exactly vanish. You will, in fact, have done theintegralexactly,andthe MonteCarlo numericalevaluationsaresuper fluous. So, backingofffromtheextremelimitingcase, to theextent thatyouareabletomake f approximatelyconstantbychangeofvariable,and totheextent thatyoucansamplea regiononlyslightlylargerthan W,youwillincreasetheaccuracyoftheMonteCarlo integral. Thistechniqueis genericallycalled reductionofvariance inthe literature. The fundamental disadvantage of simple Monte Carlo integration is that its accuracy increases only as the square root of N, the number of sampled points. If your accuracy requirements are modest, or if your computer budget is large, then the technique is highly recommended as one of great generality. In the next two sectionswe will see thattherearetechniquesavailablefor “breakingthesquareroot ofNbarrier”and achieving, at least in some cases, higher accuracy with fewer function evaluations. CITED REFERENCES AND FURTHER READING: Hammersley, J.M., and Handscomb, D.C. 1964, Monte Carlo Methods (London: Methuen). Shreider, Yu. A. (ed.) 1966, The Monte Carlo Method (Oxford: Pergamon). Sobol’, I.M. 1974, The Monte Carlo Method (Chicago: University of Chicago Press). Kalos, M.H., and Whitlock, P.A. 1986, Monte Carlo Methods (New York: Wiley). 7.7 Quasi- (that is, Sub-) Random Sequences We have just seen that choosing Npoints uniformly randomly in an n- dimensional space leads to an error term in Monte Carlo integration that decreases as1/√ N. In essence,each newpointsampledaddslinearlyto anaccumulatedsum that will become the function average, and also linearly to an accumulated sum of squares that will become the variance (equation 7.6.2). The estimated error comes from the square root of this variance, hence the power N−1/2. Just because this square root convergenceis familiar does not, however, mean that it is inevitable. A simple counterexample is to choose sample points that lie on a Cartesian grid, and to sample each grid point exactly once (in whateverorder).The MonteCarlo methodthus becomesa deterministic quadraturescheme —albeit a simple one —whose fractionalerrordecreasesat least as fast as N −1(evenfaster if the function goes to zero smoothly at the boundaries of the sampled region, or is periodic in the region).