f7-8

306 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).dominated by oneof the design parameters, that parameter will be found with this sampling technique. On the other hand, if there is an important interaction amongdifferentdesignparameters,thentheLatin hypercubegivesnoparticularadvantage. Use with care. CITED REFERENCES AND FURTHER READING: Halton, J.H. 1960, Numerische Mathematik , vol. 2, pp. 84–90. [1] Bratley P., and Fox, B.L. 1988, ACM Transactions on Mathematical Software , vol. 14, pp. 88– 100. [2] Lambert, J.P. 1988, in Numerical Mathematics – Singapore 1988 , ISNM vol. 86, R.P. Agarwal, Y.M. Chow, and S.J. Wilson, eds. (Basel: Birkha¨ user), pp. 273–284. Niederreiter, H. 1988, in Numerical Integration III , ISNM vol. 85, H. Brass and G. H¨ ammerlin, eds. (Basel: Birkha¨ user), pp. 157–171. Sobol’, I.M. 1967, USSR Computational Mathematics and Mathematical Physics , vol. 7, no. 4, pp. 86–112. [3] Antonov, I.A., and Saleev, V.M 1979, USSR Computational Mathematics and Mathematical Physics, vol. 19, no. 1, pp. 252–256. [4] Dunn, O.J., andClark, V.A. 1974, AppliedStatistics: Analysis of Variance andRegression (New York, Wiley) [discusses Latin Square]. 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). Wenow returntoitinaslightlymore formalway. Suppose that anintegrand fcanbe written astheproduct ofafunction hthatisalmostconstant timesanother, positive,function g. Then its integral over a multidimensional volume Vis/integraldisplay fd V =/integraldisplay (f/g )gdV =/integraldisplay hgd V (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 gdVa perfect differential. We then proceeded to use the basic theorem of Monte Carlo integration, equation (7.6.1). A more generalinterpretation of equation (7.8.1) is that we can integrate fby instead sampling h— not, however, with uniform probability density dV, but rather with nonuniform density gdV.I n this second interpretation, the firstinterpretation follows as the special case, where the means of generating the nonuniform sampling of gdVis 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 x iare chosen within the volume Vwith a probability density psatisfying /integraldisplay pd V =1 ( 7.8.2 ) 7.8AdaptiveandRecursiveMonteCarloMethods 307Sample 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).Thegeneralized fundamental theoremisthattheintegralofany function fisestimated,using Nsample points x1,...,x N,b y I≡/integraldisplay fd V =/integraldisplayf ppdV≈/angbracketleftbiggf p/angbracketrightbigg ±/radicalBigg /angbracketleftf2/p2/angbracketright−/angbracketleftf/p/angbracketright2 N(7.8.3 ) where angle brackets denote arithmetic means over the Npoints, exactly as in equation (7.6.2). As in equation (7.6.1), the “plus-or-minus” term is a one standard deviation errorestimate. Notice that equation (7.6.1) is in fact the special case of equation (7.8.3), withp=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/pas 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 varianceper sample point. Both angle brackets are themselves Monte Carlo estimators of integrals,so we can write S≡/angbracketleftbiggf 2 p2/angbracketrightbigg −/angbracketleftbiggf p/angbracketrightbigg2 ≈/integraldisplayf2 p2pdV−/bracketleftbigg/integraldisplayf ppdV/bracketrightbigg2 =/integraldisplayf2 pdV−/bracketleftbigg/integraldisplay fd V/bracketrightbigg2 (7.8.4 ) Wenowfindtheoptimal psubjecttotheconstraintequation(7.8.2)bythefunctionalvariation 0=δ δp/parenleftBigg/integraldisplayf2 pdV−/bracketleftbigg/integraldisplay fd V/bracketrightbigg2 +λ/integraldisplay pd V/parenrightBigg (7.8.5 ) withλa Lagrange multiplier. Note that the middle term does not depend on p. The variation (which comes inside the integrals) gives 0=−f2/p2+λor p=|f|√ λ=|f|/integraltext |f|dV(7.8.6 ) where λhas been chosen to enforce the constraint (7.8.2). Iffhas one sign in the region of integration, then we get the obvious result that the optimal choice of p— if one can figure out a practical way of effecting the sampling — is that it be proportional to |f|. Then the variance is reduced to zero. Not so obvious, but seen to be true, is the fact that p∝|f|is optimal even if ftakes on both signs. In that case the variance per sample point (from equations 7.8.4 and 7.8.6) is S=Soptimal =/parenleftbigg/integraldisplay |f|dV/parenrightbigg2 −/parenleftbigg/integraldisplay fd V/parenrightbigg2 (7.8.7 ) One curiosity is that one can add a constant to the integrand to make it all of one sign, since this changes the integral by a known amount, constant ×V. Then, the optimal choice ofpalways gives zero variance, that is, a perfectly accurate integral! The resolution of this seeming paradox (already mentioned at the end of §7.6) is that perfect knowledge of p in equation (7.8.6) requires perfect knowledge of/integraltext |f|dV, which is tantamount to already knowing the integral you are trying to compute! If your function ftakes on a known constant value in most of the volume V,i ti s certainly a good idea to add a constant so as to make that value zero. Having done that, theaccuracy attainable by importance sampling depends in practice not on how small equation(7.8.7) is, but rather on how small is equation (7.8.4) for an implementable p, likely only a crude approximation to the ideal. 308 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).Stratified Sampling The idea of stratified sampling is quite different from importance sampling. Let us expand our notation slightly and let /angbracketleft/angbracketleftf/angbracketright/angbracketrightdenote the true average of the function fover the volume V(namely the integral divided by V), while /angbracketleftf/angbracketrightdenotes as before the simplest (uniformly sampled) Monte Carlo estimator of that average: /angbracketleft/angbracketleftf/angbracketright/angbracketright ≡1 V/integraldisplay fd V /angbracketleftf/angbracketright≡1 N/summationdisplay if(xi)( 7.8.8 ) The variance of the estimator, Var (/angbracketleftf/angbracketright), which measures the square of the error of the Monte Carlo integration, is asymptotically related to the variance of the function, Var (f)≡ /angbracketleft/angbracketleftf2/angbracketright/angbracketright − /angbracketleft/angbracketleft f/angbracketright/angbracketright2, by the relation Var (/angbracketleftf/angbracketright)=Var (f) N(7.8.9 ) (compare equation 7.6.1). Suppose we divide the volume Vinto two equal, disjoint subvolumes, denoted aandb, and sample N/2points in each subvolume. Then another estimator for /angbracketleft/angbracketleftf/angbracketright/angbracketright, different from equation (7.8.8), which we denote /angbracketleftf/angbracketright/prime,i s /angbracketleftf/angbracketright/prime≡1 2/parenleftbig /angbracketleftf/angbracketrighta+/angbracketleftf/angbracketrightb/parenrightbig (7.8.10 ) in other words, the mean of the sample averages in the two half-regions. The variance of estimator (7.8.10) is given by Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig =1 4/bracketleftbig Var/parenleftbig /angbracketleftf/angbracketrighta/parenrightbig +Var/parenleftbig /angbracketleftf/angbracketrightb/parenrightbig/bracketrightbig =1 4/bracketleftbiggVara(f) N/2+Varb(f) N/2/bracketrightbigg =1 2N[Vara(f)+Varb(f)](7.8.11 ) Here Var a(f)denotes the variance of fin subregion a, that is, /angbracketleft/angbracketleftf2/angbracketright/angbracketrighta−/angbracketleft /angbracketleftf/angbracketright/angbracketright2 a, and correspondingly for b. From the definitions already given, it is not difficult to prove the relation Var (f)=1 2[Vara(f)+Varb(f)] +1 4(/angbracketleft/angbracketleftf/angbracketright/angbracketrighta−/angbracketleft /angbracketleftf/angbracketright/angbracketrightb)2(7.8.12 ) (In physics, this formula for combining second moments is the “parallel axis theorem.”) Comparing equations (7.8.9), (7.8.11), and (7.8.12), one sees that the stratified (into twosubvolumes) sampling gives a variance that is never larger than the simple Monte Carlo case— and smaller whenever the means of the stratifiedsamples, /angbracketleft/angbracketleftf/angbracketright/angbracketright aand/angbracketleft/angbracketleftf/angbracketright/angbracketrightb, are different. We have not yet exploited the possibility of sampling the two subvolumes with different numbers of points, say Nain subregion aandNb≡N−Nain subregion b. Let us do so now. Then the variance of the estimator is Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig =1 4/bracketleftbiggVara(f) Na+Varb(f) N−Na/bracketrightbigg (7.8.13 ) which is minimized (one can easily verify) when Na N=σa σa+σb(7.8.14 ) Here we have adopted the shorthand notation σa≡[Vara(f)]1/2, and correspondingly for b. IfNasatisfies equation (7.8.14), then equation (7.8.13) reduces to Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig =(σa+σb)2 4N(7.8.15 ) 7.8AdaptiveandRecursiveMonteCarloMethods 309Sample 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).Equation(7.8.15)reducestoequation(7.8.9)ifVar (f)=Vara(f)=Varb(f),inwhichcase stratifying the sample makes no difference. A standard way to generalize the above result is to consider the volume Vdivided into morethantwoequalsubregions. Onecanreadilyobtaintheresultthattheoptimalallocationofsamplepointsamong theregions istohavethenumberofpointsineach region jproportional toσ j(that is, the square root of the variance of the function fin that subregion). In spaces of high dimensionality (say d>∼4) this is not in practice very useful, however. Dividing a volume into Ksegments along each dimension implies Kdsubvolumes, typically much too large a number when one contemplates estimating all the corresponding σj’s. Mixed Strategies Importance sampling and stratified sampling seem, at first sight, inconsistent with each other. The former concentrates sample points where the magnitude of the integrand |f|is largest, that latter where the variance of fis largest. How can both be right? The answer is that (like so much else in life) it all depends on what you know and how well you know it. Importance sampling depends on already knowing some approximation toyour integral, so that you are able to generate random points x iwith the desired probability density p. To the extent that your pis not ideal, you are left with an error that decreases only as N−1/2. Things are particularly bad if your pis far from ideal in a region where the integrand fis changing rapidly, since then the sampled function h=f/pwill have a large variance. Importance samplingworksbysmoothing thevaluesofthesampledfunction h,and is effective only to the extent that you succeed in this. Stratified sampling, by contrast, does not necessarily require that you know anything about f. Stratifiedsampling works by smoothing out the fluctuations of the numberof points in subregions, not by smoothing the values of the points. The simplest stratified strategy,dividing VintoNequal subregions and choosing one point randomly in each subregion, already gives a method whose error decreases asymptotically as N −1, much faster than N−1/2. (Notethatquasi-random numbers, §7.7,areanother wayofsmoothingfluctuationsin the density of points, giving nearly as good a result as the “blind” stratification strategy.) However, “asymptotically” is an important caveat: For example, if the integrand is negligible in all but a single subregion, then the resulting one-sample integration is all butuseless. Information, even very crude, allowing importance sampling to put many points inthe active subregion would be much better than blind stratified sampling. Stratified sampling really comes into its own if you have some way of estimating the variances, so thatyou can put unequal numbers ofpoints in differentsubregions, according to(7.8.14) or its generalizations, andif you can find a way of dividing a region into a practical number of subregions (notably notK dwith large dimension d), while yet significantly reducing the variance of the function in each subregion compared to its variance in the fullvolume. Doing this requires a lot of knowledge about f, though different knowledge from what is required for importance sampling. In practice, importance sampling and stratifiedsampling are not incompatible. In many, if not most, cases of interest, the integrand fis small everywhere in Vexcept for a small fractional volume of “active regions.” In these regions the magnitude of |f|and the standard deviation σ=[Var (f)] 1/2are comparable in size, so both techniques will give about the same concentration of points. In more sophisticated implementations, it is also possible to“nest” the two techniques, so that (e.g.) importance sampling on a crude grid is followedby stratification within each grid cell. Adaptive Monte Carlo: VEGAS The VEGAS algorithm, invented by Peter Lepage [1,2], is widely used for multidimen- sional integrals that occur in elementary particle physics. VEGAS is primarily based onimportance sampling, but it also does some stratified sampling if the dimension dis small enough toavoid K dexplosion (specifically,if (K/2)d<N / 2, with Nthenumber ofsample 310 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).points). The basic technique for importance sampling in VEGAS is to construct, adaptively, a multidimensional weight function gthat isseparable , p∝g(x ,y,z,... )=gx(x)gy(y)gz(z)... (7.8.16 ) Such a function avoids the Kdexplosion in two ways: (i) It can be stored in the computer asdseparate one-dimensional functions, each defined by Ktabulated values, say — so that K×dreplaces Kd. (ii)Itcan be sampled asa probability density by consecutively sampling thedone-dimensional functions to obtain coordinate vector components (x ,y,z,... ). The optimal separable weight function can be shown to be [1] gx(x)∝/bracketleftbigg/integraldisplay dy/integraldisplay dz . . .f2(x ,y,z,... ) gy(y)gz(z).../bracketrightbigg1/2 (7.8.17 ) (and correspondingly for y,z,...). Notice that this reduces to g∝|f|(7.8.6) in one dimension. Equation (7.8.17) immediately suggests VEGAS’ adaptive strategy: Given aset of g-functions (initially all constant, say), one samples the function f, accumulating not only the overall estimator of the integral, but also the Kdestimators ( Ksubdivisions of the independent variable in each of ddimensions) of the right-hand side of equation (7.8.17). These then determine improved gfunctions for the next iteration. When the integrand fis concentrated in one, or at most a few, regions in d-space, then the weight function g’s quickly become large at coordinate values that are the projections of these regions onto the coordinate axes. The accuracy of the Monte Carlo integration is thenenormously enhanced over what simple Monte Carlo would give. The weakness of VEGAS is the obvious one: To the extent that the projection of the function fonto individual coordinate directions is uniform, VEGAS gives no concentration of sample points in those dimensions. The worst case for VEGAS, e.g., is an integrand thatis concentrated close to a body diagonal line, e.g., one from (0,0,0,... )to(1,1,1,... ). Since this geometry is completely nonseparable, VEGAS can give no advantage at all. Moregenerally, VEGAS may not do well when the integrand is concentrated in one-dimensional(or higher) curved trajectories (or hypersurfaces), unless these happen to be oriented closeto the coordinate directions. The routine vegasthat follows is essentially Lepage’s standard version, minimally modified to conform to our conventions. (We thank Lepage for permission to reproduce theprogram here.) For consistency with other versions of the VEGAS algorithm in circulation,we have preserved original variable names. The parameter NDMXis what we have called K, themaximumnumberofincrementsalongeachaxis; MXDIMisthemaximumvalueof d;some other parameters are explained in the comments. The vegasroutine performs m=itmxstatistically independent evaluations of the desired integral,each with N=ncallfunction evaluations. Whilestatistically independent, these iterations do assist each other, since each one is used to refine the sampling grid forthe next one. The results of all iterations are combined into a single best answer, and itsestimated error, by the relations I best =m/summationdisplay i=1Ii σ2 i/slashBiggm/summationdisplay i=11 σ2 iσbest =/parenleftBiggm/summationdisplay i=11 σ2 i/parenrightBigg−1/2 (7.8.18 ) Also returned is the quantity χ2/m≡1 m−1m/summationdisplay i=1(Ii−Ibest )2 σ2 i(7.8.19 ) If this is significantly larger than 1, then the results of the iterations are statistically inconsistent, and the answers are suspect. The input flag initcan be used to advantage. One might have a call with init=0, ncall=1000 ,itmx=5immediatelyfollowedbyacallwith init=1,ncall=100000 ,itmx=1. Theeffectwouldbetodevelopasamplinggridover5iterationsofasmallnumberofsamples,then to do a single high accuracy integration on the optimized grid. 7.8AdaptiveandRecursiveMonteCarloMethods 311Sample 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).Note that the user-supplied integrand function, fxn, has an argument wgtin addition to the expected evaluation point x. In most applications you ignore wgtinside the function. Occasionally,however,youmaywanttointegratesomeadditionalfunctionorfunctionsalongwith the principal function f. The integral of any such function gcan be estimated by I g=/summationdisplay iwig(x)( 7.8.20 ) where the wi’s andx’s are the arguments wgtandx, respectively. It is straightforward to accumulate this sum inside your function fxn, and to pass the answer back to your main program via a common block. Of course, g(x)had better resemble the principal function f to some degree, since the sampling will be optimized for f. SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn, * tgral,sd,chi2a) INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM REAL tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.e-30)EXTERNAL fxn C USES fxn,ran2,rebin Performs Monte Carlo integration of a user-supplied ndim -dimensional function fxn over a rectangular volume specified by region ,a 2×ndim vector consisting of ndim “lower left” coordinates of the region followed by ndim “upper right” coordinates. The integration consists of itmx iterations, each with approximately ncall calls to the function. After each iteration the grid is refined; more than 5 or 10 iterations are rarely useful. The input flag init signals whether this call is a new start, or a subsequent call for additional iterations (see comments below). The input flag nprn (normally 0) controls the amount of diagnostic output. Returned answers are tgral (the best estimate of the integral), sd(its standard deviation), and chi2a (χ2per degree of freedom, an indicator of whether consistent results are being obtained). See text for further details. INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM) REAL calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo, * d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),* r(NDMX),x(MXDIM),xi(NDMX,MXDIM),xin(NDMX),ran2 DOUBLE PRECISION schi,si,swgt COMMON /ranno/ idum Means for random number initialization. SAVE Best make everything static, allowing restarts. if(init.le.0)then Normal entry. Enter here on a cold start. mds=1 Change to mds=0 to disable stratified sampling, i.e., use im- portance sampling only. ndo=1 do 11j=1,ndim xi(1,j)=1. enddo 11 endifif (init.le.1)then Enter here to inherit the grid from a previous call, but not its answers. si=0.d0 swgt=0.d0schi=0.d0 endif if (init.le.2)then Enter here to inherit the previous grid and its answers. nd=NDMXng=1 if(mds.ne.0)then Set up for stratification. ng=(ncall/2.+0.25)**(1./ndim)mds=1 if((2*ng-NDMX).ge.0)then mds=-1npg=ng/NDMX+1nd=ng/npg ng=npg*nd endif endif k=ng**ndim 312 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).npg=max(ncall/k,2) calls=float(npg)*float(k) dxg=1./ng dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.)xnd=nd dxg=dxg*xnd xjac=1./callsdo 12j=1,ndim dx(j)=region(j+ndim)-region(j) xjac=xjac*dx(j) enddo 12 if(nd.ne.ndo)then Do binning if necessary. do13i=1,max(nd,ndo) r(i)=1. enddo 13 do14j=1,ndim call rebin(ndo/xnd,nd,r,xin,xi(1,j)) enddo 14 ndo=nd endif if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn, * ALPH,mds,nd,(j,region(j),j,region(j+ndim),j=1,ndim) endif do28it=1,itmx Main iteration loop. Can enter here ( init≥3)t od oa na d d i t i o n a l itmx iterations with all other parameters unchanged. ti=0. tsi=0. do16j=1,ndim kg(j)=1 do15i=1,nd d(i,j)=0.di(i,j)=0. enddo 15 enddo 16 10 continue fb=0.f2b=0. do 19k=1,npg wgt=xjacdo 17j=1,ndim xn=(kg(j)-ran2(idum))*dxg+1. ia(j)=max(min(int(xn),NDMX),1)if(ia(j).gt.1)then xo=xi(ia(j),j)-xi(ia(j)-1,j) rc=xi(ia(j)-1,j)+(xn-ia(j))*xo else xo=xi(ia(j),j) rc=(xn-ia(j))*xo endifx(j)=region(j)+rc*dx(j)wgt=wgt*xo*xnd enddo 17 f=wgt*fxn(x,wgt) f2=f*f fb=fb+f f2b=f2b+f2do 18j=1,ndim di(ia(j),j)=di(ia(j),j)+f if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2 enddo 18 enddo 19 f2b=sqrt(f2b*npg) f2b=(f2b-fb)*(f2b+fb) 7.8AdaptiveandRecursiveMonteCarloMethods 313Sample 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 (f2b.le.0.) f2b=TINY ti=ti+fb tsi=tsi+f2b if(mds.lt.0)then Use stratified sampling. do21j=1,ndim d(ia(j),j)=d(ia(j),j)+f2b enddo 21 endif do22k=ndim,1,-1 kg(k)=mod(kg(k),ng)+1 if(kg(k).ne.1) goto 10 enddo 22 tsi=tsi*dv2g Compute final results for this iteration. wgt=1./tsi si=si+dble(wgt)*dble(ti)schi=schi+dble(wgt)*dble(ti)**2 swgt=swgt+dble(wgt) tgral=si/swgtchi2a=max((schi-si*tgral)/(it-.99d0),0.d0)sd=sqrt(1./swgt) tsi=sqrt(tsi) if(nprn.ge.0)then write(*,201) it,ti,tsi,tgral,sd,chi2a if(nprn.ne.0)then do 23j=1,ndim write(*,202) j,(xi(i,j),di(i,j), * i=1+nprn/2,nd,nprn) enddo 23 endif endif do25j=1,ndim Refine the grid. Consult references to understand the subtlety of this procedure. The refinement is damped, to avoidrapid, destabilizing changes, and also compressed in rangeby the exponent ALPH .xo=d(1,j) xn=d(2,j)d(1,j)=(xo+xn)/2. dt(j)=d(1,j) do 24i=2,nd-1 rc=xo+xnxo=xn xn=d(i+1,j) d(i,j)=(rc+xn)/3.dt(j)=dt(j)+d(i,j) enddo 24 d(nd,j)=(xo+xn)/2. dt(j)=dt(j)+d(nd,j) enddo 25 do27j=1,ndim rc=0.do 26i=1,nd if(d(i,j).lt.TINY) d(i,j)=TINY r(i)=((1.-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPHrc=rc+r(i) enddo 26 call rebin(rc/xnd,nd,r,xin,xi(1,j)) enddo 27 enddo 28 return 200 FORMAT(/’ input parameters for vegas: ndim=’,i3,’ ncall=’,f8.0* /28x,’ it=’,i5,’ itmx=’,i5* /28x,’ nprn=’,i3,’ alph=’,f5.2/28x,’ mds=’,i3,’ nd=’,i4 * /(30x,’xl(’,i2,’)= ’,g11.4,’ xu(’,i2,’)= ’,g11.4)) 201 FORMAT(/’ iteration no.’,I3,’: ’,’integral =’,g14.7,’+/- ’,g9.2* /’ all iterations: integral =’,g14.7,’+/- ’,g9.2,* ’ chi**2/it’’n =’,g9.2) 202 FORMAT(/’ data for axis ’,I2/’ X delta i ’, 314 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).* ’ x delta i ’,’ x delta i ’, * /(1x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4)) END SUBROUTINE rebin(rc,nd,r,xin,xi) INTEGER ndREAL rc,r(*),xi(*),xin(*) Utility routine used by vegas , to rebin a vector of densities xiinto new bins defined by a vector r. INTEGER i,kREAL dr,xn,xo k=0 xo=0.dr=0. do 11i=1,nd-1 1 if(rc.gt.dr)then k=k+1dr=dr+r(k) goto 1 endifif(k.gt.1) xo=xi(k-1) xn=xi(k) dr=dr-rcxin(i)=xn-(xn-xo)*dr/r(k) enddo 11 do12i=1,nd-1 xi(i)=xin(i) enddo 12 xi(nd)=1. return END Recursive Stratified Sampling The problem with stratified sampling, we have seen, is that it may not avoid the Kd explosion inherent in the obvious, Cartesian, tessellation of a d-dimensional volume. A technique called recursive stratified sampling [3]attempts to do this by successive bisections of a volume, not along all ddimensions, but rather along only one dimension at a time. The starting points are equations (7.8.10) and (7.8.13), applied to bisections of successivelysmaller subregions. Suppose that we have a quota of Nevaluations of the function f, and want to evaluate /angbracketleftf/angbracketright/primein the rectangular parallelepiped region R=(xa,xb). (We denote such a region by the two coordinate vectors of its diagonally opposite corners.) First, we allocate a fraction pof Ntowards exploring the variance of finR: We sample pNfunction values uniformly in Rand accumulate the sums that will give the ddifferent pairs of variances corresponding to theddifferent coordinate directions along which Rcan be bisected. In other words, in pN samples, we estimate Var (f)in each of the regions resulting from a possible bisection of R, Rai≡(xa,xb−1 2ei·(xb−xa)ei) Rbi≡(xa+1 2ei·(xb−xa)ei,xb)(7.8.21 ) Hereeiis the unit vector in the ith coordinate direction, i=1,2,...,d. Second, we inspect the variances to find the most favorable dimension ito bisect. By equation (7.8.15), we could, for example, choose that ifor which the sum of the square roots ofthevariance estimatorsinregions RaiandRbiisminimized. (Actually,aswewillexplain, we do something slightly different.) 7.8AdaptiveandRecursiveMonteCarloMethods 315Sample 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).Third, we allocate the remaining (1−p)Nfunction evaluations between the regions RaiandRbi. Ifweused equation (7.8.15) tochoose i,weshould dothisallocation according to equation (7.8.14). Wenowhavetwoparallelepipedseachwithitsownallocationoffunctionevaluationsfor estimatingthemean of f. Our“RSS”algorithm nowshows itselftobe recursive: Toevaluate the mean in each region, we go back to the sentence beginning “First,...” in the paragraphabove equation (7.8.21). (Of course, when the allocation of points to a region falls belowsome number, we resort to simple Monte Carlo rather than continue with the recursion.) Finally, we combine the means, and also estimated variances of the two subvolumes, using equation (7.8.10) and the first line of equation (7.8.11). This completes the RSS algorithm in its simplest form. Before we describe some additionaltricksunderthegeneralrubricof“implementationdetails,”weneedtoreturnbrieflyto equations (7.8.13)–(7.8.15) and derive the equations that we actually use instead of these.The right-hand side of equation (7.8.13) applies the familiar scaling law of equation (7.8.9)twice, once to aand again to b. This would be correct if the estimates /angbracketleftf/angbracketright aand/angbracketleftf/angbracketrightbwere each made by simple Monte Carlo, with uniformly random sample points. However, the two estimatesofthemeanareinfactmaderecursively. Thus,thereisnoreasontoexpectequation(7.8.9) to hold. Rather, we might substitute for equation (7.8.13) the relation, Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig =1 4/bracketleftbiggVara(f) Nαa+Varb(f) (N−Na)α/bracketrightbigg (7.8.22 ) where αis an unknown constant ≥1(the case of equality corresponding to simple Monte Carlo). In that case, a short calculation shows that Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig is minimized when Na N=Vara(f)1/(1+α) Vara(f)1/(1+α)+Varb(f)1/(1+α)(7.8.23 ) and that its minimum value is Var/parenleftbig /angbracketleftf/angbracketright/prime/parenrightbig ∝/bracketleftBig Vara(f)1/(1+α)+Varb(f)1/(1+α)/bracketrightBig1+α (7.8.24 ) Equations (7.8.22)–(7.8.24) reduce to equations (7.8.13)–(7.8.15) when α=1. Numerical experiments to find a self-consistent value for αfind that α≈2. That is, when equation (7.8.23)with α=2isusedrecursivelytoallocatesampleopportunities,theobservedvariance of the RSS algorithm goes approximately as N−2, while any other value of αin equation (7.8.23) gives a poorer fall-off. (The sensitivity to αis, however, not very great; it is not known whether α=2is an analytically justifiable result, or only a useful heuristic.) Turn now to the routine, miser, which implements the RSS method. A bit of FORTRAN wizardry is its implementation of the required recursion. This is done by dimensioning anarray stack,andashorter“stackframe” stf;thelatterhascomponents thatareequivalenced to variables that need to be preserved during the recursion, including a flag indicating whereprogram control should return. A recursive call then consists of copying the stack frameonto the stack, incrementing the stack pointer jstack, and transferring control. A recursive return analogously pops the stack and transfers control to the saved location. Stack growthinmiseris only logarithmic in N, since at each bifurcation one of the subvolumes can be processed immediately. Theprincipaldifferencebetween miser’simplementationandthealgorithmasdescribed thusfarliesinhowthevariancesontheright-handsideofequation(7.8.23)areestimated. Wefindempiricallythatitissomewhatmorerobusttousethesquareofthedifferenceofmaximumandminimumsampledfunctionvalues,insteadofthegenuinesecondmomentofthesamples.Thisestimatorisofcourseincreasinglybiasedwithincreasingsamplesize;however,equation(7.8.23)usesitonlytocomparetwosubvolumes( aandb)havingapproximatelyequalnumbers of samples. The “max minus min” estimator proves its worth when the preliminary samplingyields only a single point, or small number of points, in active regions of the integrand. Inmany realistic cases, these are indicators of nearby regions of even greater importance, and itis useful to let them attract the greater sampling weight that “max minus min” provides. Asecondmodificationembodiedinthecodeistheintroductionofa“ditheringparameter,” dith,whosenonzero valuecausessubvolumestobedividednotexactlydownthemiddle,but 316 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).rather into fractions 0.5±dith, with the sign of the ±randomly chosen by a built-in random number routine. Normally dithcan be set to zero. However, there is a large advantage in taking dithto be nonzero if some special symmetry of the integrand puts the active region exactly at the midpoint of the region, or at the center of some power-of-two submultiple ofthe region. One wants to avoid the extreme case of the active region being evenly dividedinto 2 dabutting corners of a d-dimensional space. A typical nonzero value of dith,o n those occasions when it is useful, might be 0.1. Of course, when the dithering parameter is nonzero, we must take the differing sizes of the subvolumes into account; the code doesthis through the variable fracl. One final feature in the code deserves mention. The RSS algorithm uses a single set of sample points to evaluate equation (7.8.23) in all ddirections. At bottom levels of the recursion, the number of sample points can be quite small. Although rare, it can happen thatin one direction all the samples are in one half of the volume; in that case, that directionis ignored as a candidate for bifurcation. Even more rare is the possibility that all of thesamples are in one half of the volume in alldirections. In this case, a random direction is chosen. If this happens too often in your application, then you should increase MNPT(see lineif (jb.eq.0) ...in the code). Note that miser, as given, returns as avean estimate of the average function value /angbracketleft/angbracketleftf/angbracketright/angbracketright,not the integralof fover theregion. Theroutine vegas,adopting the other convention, returns as tgraltheintegral. Thetwo conventions are ofcourse triviallyrelated,by equation (7.8.8), since the volume Vof the rectangular region is known. SUBROUTINE miser(func,region,ndim,npts,dith,ave,var) INTEGER ndim,npts,MNPT,MNBS,MAXD,NSTACKREAL ave,dith,var,region(2*ndim),func,TINY,BIG,PFAC PARAMETER (MNPT=15,MNBS=4*MNPT,MAXD=10,TINY=1.e-30,BIG=1.e30, * NSTACK=1000,PFAC=0.1) EXTERNAL func C USES func,ranpt Monte Carlo samples a user-supplied ndim -dimensional function func in a rectangular volume specified by region ,a 2×ndim vector consisting of ndim “lower-left” coordinates of the region followed by ndim “upper-right” coordinates. The function is sampled a total ofnpts times, at locations determined by the method of recursive stratified sampling. The mean value of the function in the region is returned as ave; an estimate of the statistical uncertainty of ave (square of standard deviation) is returned as var. The input parameter dith should normally be set to zero, but can be set to (e.g.) 0.1 if func ’s active region falls on the boundary of a power-of-two subdivision of region . Parameters: PFAC is the fraction of remaining function evaluations used at each stage to explore the variance of func .A tl e a s t MNPT function evaluations are performed in any terminal subregion; a subregion is further bisected only if at least MNBS function evaluations are available. MAXD is the largest value of ndim .NSTACK is the total size of the stack. INTEGER iran,j,jb,jstack,n,naddr,np,npre,nptl,nptr,npttREAL avel,fracl,fval,rgl,rgm,rgr,s,sigl,siglb,sigr,sigrb,sum, * sumb,summ,summ2,varl,fmaxl(MAXD),fmaxr(MAXD),fminl(MAXD), * fminr(MAXD),pt(MAXD),rmid(MAXD),stack(NSTACK),stf(9) EQUIVALENCE (stf(1),avel),(stf(2),varl),(stf(3),jb), * (stf(4),nptr),(stf(5),naddr),(stf(6),rgl),(stf(7),rgm), * (stf(8),rgr),(stf(9),fracl) SAVE iranDATA iran /0/ jstack=0 nptt=npts 1 continue if (nptt.lt.MNBS) then Too few points to bisect; do straight Monte Carlo. np=abs(nptt)summ=0.summ2=0. do 11n=1,np call ranpt(pt,region,ndim)fval=func(pt) summ=summ+fval 7.8AdaptiveandRecursiveMonteCarloMethods 317Sample 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).summ2=summ2+fval**2 enddo 11 ave=summ/npvar=max(TINY,(summ2-summ**2/np)/np**2) else Do the preliminary (uniform) sampling. npre=max(int(nptt*PFAC),MNPT) do 12j=1,ndim Initialize the left and right bounds for each dimension. iran=mod(iran*2661+36979,175000)s=sign(dith,float(iran-87500)) rmid(j)=(0.5+s)*region(j)+(0.5-s)*region(j+ndim) fminl(j)=BIGfminr(j)=BIGfmaxl(j)=-BIG fmaxr(j)=-BIG enddo 12 do14n=1,npre Loop over the points in the sample. call ranpt(pt,region,ndim) fval=func(pt)do 13j=1,ndim Find the left and right bounds for each dimension. if(pt(j).le.rmid(j))then fminl(j)=min(fminl(j),fval) fmaxl(j)=max(fmaxl(j),fval) else fminr(j)=min(fminr(j),fval) fmaxr(j)=max(fmaxr(j),fval) endif enddo 13 enddo 14 sumb=BIG Choose which dimension jbto bisect. jb=0 siglb=1. sigrb=1.do 15j=1,ndim if(fmaxl(j).gt.fminl(j).and.fmaxr(j).gt.fminr(j))then sigl=max(TINY,(fmaxl(j)-fminl(j))**(2./3.)) sigr=max(TINY,(fmaxr(j)-fminr(j))**(2./3.))sum=sigl+sigr Equation (7.8.24), see text. if (sum.le.sumb) then sumb=sum jb=jsiglb=sigl sigrb=sigr endif endif enddo 15 if (jb.eq.0) jb=1+(ndim*iran)/175000 MNPT may be too small. rgl=region(jb) Apportion the remaining points between left and right. rgm=rmid(jb) rgr=region(jb+ndim) fracl=abs((rgm-rgl)/(rgr-rgl))nptl=MNPT+(nptt-npre-2*MNPT) * *fracl*siglb/(fracl*siglb+(1.-fracl)*sigrb) Equation (7.8.23). nptr=nptt-npre-nptl region(jb+ndim)=rgm Set region to left. naddr=1 Push the stack. do 16j=1,9 stack(jstack+j)=stf(j) enddo 16 jstack=jstack+9 nptt=nptl goto 1 Dispatch recursive call; will return back here eventually. 10 continue avel=ave Save left estimates on stack variable. varl=var 318 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).region(jb)=rgm Set region to right. region(jb+ndim)=rgr naddr=2 Push the stack. do17j=1,9 stack(jstack+j)=stf(j) enddo 17 jstack=jstack+9 nptt=nptrgoto 1 Dispatch recursive call; will return back here eventually. 20 continue region(jb)=rgl Restore region to original value (so that we don’t need to include it on the stack). ave=fracl*avel+(1.-fracl)*ave var=fracl**2*varl+(1.-fracl)**2*var Combine left and right regions by equa- tion (7.8.11) (1st line). endif if (jstack.ne.0) then Pop the stack. jstack=jstack-9 do 18j=1,9 stf(j)=stack(jstack+j) enddo 18 goto (10,20),naddr pause ’miser: never get here’ endifreturn END Themiserroutinecallsashortsubroutine ranpttogetarandompointwithinaspecified d-dimensional region. The following version of ranptmakes consecutive calls to a uniform random number generator and does the obvious scaling. One can easily modify ranptto generate its points via the quasi-random routine sobseq(§7.7). We find that miserwith sobseqcan be considerably more accurate than miserwith uniform random deviates. Since the use of RSS and the use of quasi-random numbers are completely separable, however, wehave not made the code given here dependent on sobseq. A similar remark might be made regarding importance sampling, which could in principle be combined with RSS.(One couldin principle combine vegasandmiser, although the programming would be intricate.) SUBROUTINE ranpt(pt,region,n) INTEGER n,idum REAL pt(n),region(2*n)COMMON /ranno/ idumSAVE /ranno/ C USES ran1 Returns a uniformly random point ptin an n-dimensional rectangular region .U s e db y miser ; calls ran1 for uniform deviates. Your main program should initialize idum , through theCOMMON block /ranno/ , to a negative seed integer. INTEGER jREAL ran1do 11j=1,n pt(j)=region(j)+(region(j+n)-region(j))*ran1(idum) enddo 11 return END CITED REFERENCES AND FURTHER READING: Hammersley, J.M. and Handscomb, D.C. 1964, Monte Carlo Methods (London: Methuen). Kalos, M.H. and Whitlock, P.A. 1986, Monte Carlo Methods (New York: Wiley). Bratley,P.,Fox,B.L.,andSchrage,E.L.1983, AGuidetoSimulation (NewYork:Springer-Verlag). Lepage, G.P. 1978, Journal of Computational Physics , vol. 27, pp. 192–203. [1] 7.8AdaptiveandRecursiveMonteCarloMethods 319Sample 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).Lepage, G.P. 1980, “VEGAS: An Adaptive Multidimensional Integration Program,” Publication CLNS-80/447, Cornell University. [2] Press, W.H., and Farrar, G.R. 1990, Computers in Physics , vol. 4, pp. 190–195. [3]