f7-7

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 numericalevaluationsaresuperfluous. 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“breakingthesquarerootofNbarrier” 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). 300 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).The trouble with a grid is that one has to decide in advance how fine it should be. One is then committed to completing all of its sample points. With a grid, it isnot convenient to “sample until” some convergence or termination criterion is met. One might ask if there is not some intermediate scheme, some way to pick sample points “at random,” yet spread out in some self-avoiding way, avoiding the chanceclustering that occurs with uniformly random points. AsimilarquestionarisesfortasksotherthanMonteCarlointegration. Wemight want to searchan n-dimensionalspace fora pointwheresome (locallycomputable) condition holds. Of course, for the task to be computationally meaningful, there had better be continuity, so that the desired condition will hold in some finite n- dimensionalneighborhood. We maynotknow a priorihowlargethatneighborhood is, however. We wantto“sample until”thedesiredpointis found,movingsmoothly to finer scales with increasing samples. Is there any way to do this that is betterthan uncorrelated, random samples? The answer to the above question is “yes.” Sequences of n-tuples that fill n-space more uniformly than uncorrelated random points are called quasi-random sequences . That term is somewhat of a misnomer, since there is nothing “random” aboutquasi-randomsequences: Theyarecleverlycraftedtobe,infact, sub-random. The sample points in a quasi-random sequence are, in a precise sense, “maximally avoiding” of each other. A conceptuallysimpleexampleis Halton’ssequence [1]. Inonedimension,the jth number Hjin the sequence is obtained by the following steps: (i) Write jas a numberinbase b,where bissomeprime. (Forexample j=1 7inbase b=3is122.) (ii) Reverse the digits and put a radix point (i.e., a decimal point base b) in front of the sequence. (In the example, we get 0.221base 3.) The result is H j. To get a sequenceof n-tuplesin n-space,youmakeeachcomponenta Halton sequencewith a different prime base b. Typically, the first nprimes are used. It is not hard to see how Halton’s sequence works: Every time the number of digits in jincreases by one place, j’s digit-reversed fraction becomes a factor of bfiner-meshed. Thus the process is one of filling in all the points on a sequence of finer and finer Cartesian grids — and in a kind of maximally spread-out order on each grid (since, e.g., the most rapidly changing digit in jcontrols the most significant digit of the fraction). Other ways of generating quasi-random sequences have been suggested by Faure, Sobol’, Niederreiter, and others. Bratley and Fox [2]provide a good review andreferences,anddiscuss a particularlyefficientvariantof the Sobol’ [3]sequence suggested by Antonov and Saleev [4]. It is this Antonov-Saleev variant whose implementation we now discuss. TheSobol’sequencegeneratesnumbersbetweenzeroandonedirectlyasbinaryfractions oflength wbits,fromasetof wspecialbinaryfractions, Vi,i=1,2,...,w,calleddirection numbers. In Sobol’s original method, the jth number Xjis generated by XORing (bitwise exclusive or) together the set of Vi’ssatisfying the criterion on i, “the ith bit of jis nonzero.” Asjincrements, in other words, different ones of the Vi’s flash in and out of Xjon different timescales. V1alternatesbetweenbeingpresentandabsentmostquickly, while Vkgoesfrom present to absent (or vice versa) only every 2k−1steps. Antonov and Saleev’s contribution was to show that instead of using the bits of the integer jto select direction numbers, one could just as well use the bits of the Gray code of j,G(j). (For a quick review of Gray codes, look at §20.2.) NowG(j)andG(j+1 )differ in exactly one bit position, namely in the position of the 7.7Quasi-(thatis, Sub-)RandomSequences 301Sample 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).. ....... ......... .... ..... .... .............. ......... ... ............. .................... ................... .... ................ 0. 2. 4. 6. 81 points 1 to 1280.2.4.6.81 ........ ......... .... ..... . ... ....... ... .... ............ ............. . . ..... .. ........... ............. ...... .... ............................................................................................................................. ........... .................. .... ......................................................................... ..................... ................... 0. 2. 4. 6. 81 points 129 to 5120.2.4.6.81 ............................................. ......................................................................................................................................................................................................................................................................................................................................................................................... ............... ............................................. 0. 2. 4. 6. 81 points 513 to 10240.2.4.6.81 ............ .............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................. ................................................................................................................. .............................................................................................................................................................................................................................................................................................................. 0. 2. 4. 6. 81 points 1 to 10240.2.4.6.81 Figure 7.7.1. First 1024 points of a two-dimensional Sobol ’sequence. The sequence is generated number-theoretically, rather than randomly, so successive points at any stage “know”how tofill in the gaps in the previously generated distribution. rightmostzerobitinthebinaryrepresentationof j(addingaleadingzeroto jifnecessary). A consequence is that the j+1st Sobol’-Antonov-Saleev number can be obtained from the jth by XORing it with a single Vi, namely with ithe position of the rightmost zero bit in j. This makes the calculation of the sequence very ef ficient, as we shall see. Figure7.7.1plotsthe first1024pointsgeneratedbyatwo-dimensionalSobol ’sequence. One sees that successive points do “know”about the gaps left previously, and keep filling them in, hierarchically. Wehavedeferredtothispointadiscussionofhowthedirectionnumbers Viaregenerated. Some nontrivial mathematics is involved in that, so we will content ourself with a cookbooksummaryonly: EachdifferentSobol ’sequence (orcomponentofan n-dimensionalsequence) is based on a different primitive polynomial over the integers modulo 2, that is, a polynomialwhose coef ficients are either 0 or 1, and which generates a maximal length shift register sequence. (Primitive polynomials modulo 2 were used in §7.4, and are further discussed in §20.3.) Suppose Pis such a polynomial, of degree q, P=x q+a1xq−1+a2xq−2+··· +aq−1x+1 ( 7.7.1 ) 302 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).Degree Primitive Polynomials Modulo 2* 10 (i.e., x+1) 21 (i.e., x2+x+1) 31, 2 (i.e., x3+x+1and x3+x2+1) 41, 4 (i.e., x4+x+1and x4+x3+1) 52, 4, 7, 11, 13, 14 61, 13, 16, 19, 22, 25 71, 4, 7, 8, 14, 19, 21, 28, 31, 32, 37, 41, 42, 50, 55, 56, 59, 62 814, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84, 97, 103, 115, 122 98, 13, 16, 22, 25, 44, 47, 52, 55, 59, 62, 67, 74, 81, 82, 87, 91, 94, 103, 104, 109, 122, 124, 137, 138, 143, 145, 152, 157, 167, 173, 176, 181, 182, 185, 191, 194, 199, 218, 220, 227, 229, 230, 234, 236, 241, 244, 253 104, 13, 19, 22, 50, 55, 64, 69, 98, 107, 115, 121, 127, 134, 140, 145, 152, 158, 161, 171, 181, 194, 199, 203, 208, 227, 242, 251, 253, 265, 266, 274, 283, 289, 295, 301, 316,319, 324, 346, 352, 361, 367, 382, 395, 398, 400, 412, 419, 422, 426, 428, 433, 446,454, 457, 472, 493, 505, 508 *Expressed as a decimal integer representing the interior bits (that is, omitting the high-order bit and the unit bit). Define a sequence of integers M iby the q-term recurrence relation, M i=2a1M i−1⊕22a2M i−2⊕···⊕ 2q−1M i−q+1aq−1⊕(2qM i−q⊕M i−q)(7.7.2 ) HerebitwiseXORisdenotedby ⊕. Thestartingvaluesforthisrecurrencearethat M1,...,M q can be arbitrary odd integers less than 2,..., 2q, respectively. Then, the direction numbers Viare given by Vi=M i/2ii=1,...,w (7.7.3 ) The accompanying table lists all primitive polynomials modulo 2 with degree q≤10. Sincethecoef ficientsareeither0or1,andsincethecoef ficientsof xqandof 1arepredictably 1,itisconvenienttodenoteapolynomialbyitsmiddlecoef ficientstakenasthebitsofabinary number (higher powers of xbeing more signi ficant bits). The table uses this convention. TurnnowtotheimplementationoftheSobol ’sequence. Successivecallstothefunction sobseq(after a preliminary initializing call) return successive points in an n-dimensional Sobol’sequence based on the firstnprimitive polynomials in the table. As given, the routine is initialized for maximum nof 6 dimensions, and for a word length wof 30 bits. These parameters can be altered by changing MAXBIT(≡w) and MAXDIM, and by adding more initializing data to the arrays ip(the primitive polynomials from the table), mdeg(their degrees), and iv(the starting values for the recurrence, equation 7.7.2). A second table, below, elucidates the initializing data in the routine. SUBROUTINE sobseq(n,x) INTEGER n,MAXBIT,MAXDIM REAL x(*) PARAMETER (MAXBIT=30,MAXDIM=6) When nis negative, internally initializes a set of MAXBITdirection numbers for each of MAXDIMdifferent Sobol’ sequences. When nis positive (but ≤MAXDIM), returns as the 7.7Quasi-(thatis,Sub-)RandomSequences 303Sample 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).Initializing Values Used in sobseq Degree Polynomial Starting Values 1 0 1(3)(5)(15) ... 2 1 11(7)(11) ... 3 1 137(5) ... 3 2 133(15) ... 4 1 11313 ... 4 4 1159 ... Parenthesized values are not freely speci fiable, but are forced by the required recurrence for this degree. vector x(1..n)thenextvaluesfrom nofthesesequences. ( nmustnotbechangedbetween initializations.) INTEGER i,im,in,ipp,j,k,l,ip(MAXDIM),iu(MAXDIM,MAXBIT), * iv(MAXBIT*MAXDIM),ix(MAXDIM),mdeg(MAXDIM) REAL fac SAVE ip,mdeg,ix,iv,in,facEQUIVALENCE (iv,iu) To allowboth 1D and 2D addressing. DATA ip /0,1,1,2,1,4/, mdeg /1,2,3,3,4,4/, ix /6*0/ DATA iv /6*1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9,156*0/ if (n.lt.0) then Initialize, don’t return a vector. do 11k=1,MAXDIM ix(k)=0 enddo 11 in=0if(iv(1).ne.1)return fac=1./2.**MAXBIT do 15k=1,MAXDIM do12j=1,mdeg(k) Stored values only require normalization. iu(k,j)=iu(k,j)*2**(MAXBIT-j) enddo 12 do14j=mdeg(k)+1,MAXBIT Use the recurrence to get other values. ipp=ip(k) i=iu(k,j-mdeg(k)) i=ieor(i,i/2**mdeg(k))do 13l=mdeg(k)-1,1,-1 if(iand(ipp,1).ne.0)i=ieor(i,iu(k,j-l)) ipp=ipp/2 enddo 13 iu(k,j)=i enddo 14 enddo 15 else Calculate the next vector in the sequence. im=in do16j=1,MAXBIT Find the rightmost zero bit. if(iand(im,1).eq.0)goto 1im=im/2 enddo 16 pause ’MAXBIT too small in sobseq’ 1 im=(j-1)*MAXDIM do17k=1,min(n,MAXDIM) XORtheappropriatedirectionnumberintoeachcom- ponent of the vector and convert to a floating number.ix(k)=ieor(ix(k),iv(im+k)) x(k)=ix(k)*fac enddo 17 in=in+1 Increment the counter. 304 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).endif return END How good is a Sobol ’sequence, anyway? For Monte Carlo integration of a smooth function in ndimensions, the answer is that the fractional error will decrease with N, the numberofsamples,as (lnN)n/N,i.e.,almostasfastas 1/N. Asanexample,letusintegrate a function that is nonzero inside a torus (doughnut) in three-dimensional space. If the majorradius of the torus is R 0, the minor radial coordinate ris defined by r=/parenleftBig [(x2+y2)1/2−R0]2+z2/parenrightBig1/2 (7.7.4 ) Let us try the function f(x, y, z )=  1+c o s/parenleftbiggπr2 a2/parenrightbigg r<r 0 0 r≥r0(7.7.5 ) which can be integrated analytically in cylindrical coordinates, giving /integraldisplay/integraldisplay/integraldisplay dx dy dz f (x, y, z )=2 π2a2R0 (7.7.6 ) With parameters R0=0.6,r 0=0.3, we did 100 successive Monte Carlo integrations of equation (7.7.4), sampling uniformly in the region −1< x ,y,z < 1, for the two cases of uncorrelated randompointsandtheSobol ’sequence generated bytheroutine sobseq. Figure 7.7.2 shows the results, plotting the r.m.s. average error of the 100 integrations as a functionof the number of points sampled. (For any singleintegration, the error of course wanders from positive to negative, or vice versa, so a logarithmic plot of fractional error is not veryinformative.) The thin, dashed curve corresponds to uncorrelated random points and shows the familiar N −1/2asymptotics. The thin, solid gray curve shows the result for the Sobol ’ sequence. The logarithmic term in the expected (lnN)3/Nis readily apparent as curvature in the curve, but the asymptotic N−1is unmistakable. Tounderstand theimportance ofFigure 7.7.2,suppose that aMonte Carlointegrationof fwith1%accuracy isdesired. TheSobol ’sequence achieves thisaccuracy inafewthousand samples, while pseudorandom sampling requires nearly 100,000 samples. The ratio wouldbe even greater for higher desired accuracies. A different, not quite so favorable, case occurs when the function being integrated has hard (discontinuous) boundaries inside the sampling region, for example the function that isone inside the torus, zero outside, f(x, y, z )=/braceleftBig1 r<r 0 0 r≥r0(7.7.7 ) where risdefinedinequation(7.7.4). Notbycoincidence, thisfunctionhasthesameanalytic integral as the function of equation (7.7.5), namely 2π2a2R0. The carefully hierarchical Sobol ’sequence is based on a set of Cartesian grids, but the boundaryofthetorushasnoparticularrelationtothosegrids. Theresultisthatitisessentiallyrandom whether sampled points in a thin layer at the surface of the torus, containing on the orderof N 2/3points,comeouttobeinside,oroutside,thetorus. Thesquarerootlaw,applied to this thin layer, gives N1/3fluctuations in the sum, or N−2/3fractional error in the Monte Carlo integral. One sees this behavior veri fied in Figure 7.7.2 by the thicker gray curve. The thickerdashedcurveinFigure7.7.2istheresultofintegratingthefunction ofequation(7.7.7)using independent random points. Whilethe advantage of the Sobol ’sequence is not quite so dramatic as in the case of a smooth function, it can nonetheless be a signi ficant factor ( ∼5) even at modest accuracies like 1%, and greater at higher accuracies. Note that we have not provided the routine sobseqwith a means of starting the sequence at a point other than the beginning, but this feature would be easy to add. Oncethe initialization of the direction numbers ivhas been done, the jth point can be obtained directly by XORing together those direction numbers corresponding to nonzero bits in theGray code of j, as described above. 7.7Quasi-(thatis,Sub-)RandomSequences 305Sample 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)..001.01.1fractional accuracy of integral number of points N100 1000 10000 105∝N−1/2 ∝N−2/3 ∝N−1 quasi-random, hard boundarypseudo-random, soft boundarypseudo-random, hard boundary quasi-random, soft boundary Figure7.7.2. Fractional accuracy ofMonteCarlo integrations asafunction ofnumberofpoints sampled, for two different integrands and two different methods of choosing random points. The quasi-randomSobol’sequence converges much more rapidly than a conventional pseudo-random sequence. Quasi- random sampling does better when the integrand is smooth ( “soft boundary ”) than when it has step discontinuities ( “hard boundary ”). The curves shown are the r.m.s. average of 100 trials. The Latin Hypercube We mightheregivepassingmentiontheunrelatedtechniqueof Latinsquare or Latinhypercube sampling,whichisusefulwhenyoumustsamplean N-dimensional spaceexceedingly sparsely, at Mpoints. For example, you may want to test the crashworthinessof cars as a simultaneousfunctionof4 differentdesign parameters, but with a budget of only three expendable cars. (The issue is not whether this is a good plan —it isn’t—but rather how to make the best of the situation!) Theideaistopartitioneachdesignparameter(dimension)into Msegments,so that the whole space is partitioned into MNcells. (You can choose the segments in eachdimensiontobeequalorunequal,accordingtotaste.) With4parametersand3 cars, for example, you end up with 3×3×3×3=8 1cells. Next, choose Mcells to contain the sample points by the following algorithm: Randomly choose one of the MNcells for the first point. Now eliminate all cells that agree with this point on anyof its parameters (that is, cross out all cells in the same row, column, etc.), leaving (M−1)Ncandidates. Randomly choose one of these, eliminate new rows and columns, and continuethe process until there is only one cell left, which then contains the final sample point. The result of this construction is that eachdesign parameter will have been tested in every one of its subranges. If the response of the system under test is 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 inde finite 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 )