c7-7 (779521), страница 3
Текст из файла (страница 3)
The curves shown are the r.m.s. average of 100 trials.sequence. The logarithmic term in the expected (ln N )3 /N is readily apparent as curvaturein the curve, but the asymptotic N −1 is unmistakable.To understand the importance of Figure 7.7.2, suppose that a Monte Carlo integration off with 1% accuracy is desired. The Sobol’ sequence achieves this accuracy in a few thousandsamples, 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 hashard (discontinuous) boundaries inside the sampling region, for example the function that isone inside the torus, zero outside,n1r < r0f (x, y, z) =(7.7.7)0r ≥ r0where r is defined in equation (7.7.4).
Not by coincidence, this function has the same analyticintegral as the function of equation (7.7.5), namely 2π2 a2 R0 .The carefully hierarchical Sobol’ sequence is based on a set of Cartesian grids, but theboundary of the torus has no particular relation to those grids. The result is that it is essentiallyrandom whether sampled points in a thin layer at the surface of the torus, containing on theorder of N 2/3 points, come out to be inside, or outside, the torus. The square root law, appliedto this thin layer, gives N 1/3 fluctuations in the sum, or N −2/3 fractional error in the MonteCarlo integral. One sees this behavior verified in Figure 7.7.2 by the thicker gray curve.
Thethicker dashed curve in Figure 7.7.2 is the result of integrating the function of equation (7.7.7)using independent random points. While the advantage of the Sobol’ sequence is not quite sodramatic as in the case of a smooth function, it can nonetheless be a significant factor (∼5)even at modest accuracies like 1%, and greater at higher accuracies.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).fractional accuracy of integral.17.7 Quasi- (that is, Sub-) Random Sequences315Note that we have not provided the routine sobseq with a means of starting thesequence at a point other than the beginning, but this feature would be easy to add. Oncethe initialization of the direction numbers iv has been done, the jth point can be obtaineddirectly by XORing together those direction numbers corresponding to nonzero bits in theGray code of j, as described above.We might here give passing mention the unrelated technique of Latin square orLatin hypercube sampling, which is useful when you must sample an N -dimensionalspace exceedingly sparsely, at M points.
For example, you may want to test thecrashworthiness of cars as a simultaneous function of 4 different design parameters,but with a budget of only three expendable cars. (The issue is not whether this is agood plan — it isn’t — but rather how to make the best of the situation!)The idea is to partition each design parameter (dimension) into M segments, sothat the whole space is partitioned into M N cells. (You can choose the segments ineach dimension to be equal or unequal, according to taste.) With 4 parameters and 3cars, for example, you end up with 3 × 3 × 3 × 3 = 81 cells.Next, choose M cells to contain the sample points by the following algorithm:Randomly choose one of the M N cells for the first point.
Now eliminate all cellsthat agree with this point on any of its parameters (that is, cross out all cells in thesame row, column, etc.), leaving (M − 1)N candidates. Randomly choose one ofthese, eliminate new rows and columns, and continue the process until there is onlyone cell left, which then contains the final sample point.The result of this construction is that each design parameter will have beentested in every one of its subranges. If the response of the system under test isdominated by one of the design parameters, that parameter will be found withthis sampling technique.
On the other hand, if there is an important interactionamong different design parameters, then the Latin hypercube gives no particularadvantage. 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üser), pp. 273–284.Niederreiter, H. 1988, in Numerical Integration III, ISNM vol. 85, H. Brass and G. Hämmerlin,eds. (Basel: Birkhaüser), 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 MathematicalPhysics, vol. 19, no. 1, pp. 252–256. [4]Dunn, O.J., and Clark, V.A. 1974, Applied Statistics: Analysis of Variance and Regression (NewYork, Wiley) [discusses Latin Square].Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The Latin Hypercube316Chapter 7.Random Numbers7.8 Adaptive and Recursive Monte CarloMethodsImportance SamplingThe 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 writtenas the product of a function h that is almost constant times another, positive, function g. Thenits integral over a multidimensional volume V isZZZf 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 toG, the indefinite integral of g.
That made gdV a perfect differential. We then proceededto 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 f by instead sampling h — not,however, with uniform probability density dV , but rather with nonuniform density gdV .
Inthis second interpretation, the first interpretation follows as the special case, where the meansof generating the nonuniform sampling of gdV is via the transformation method, using theindefinite integral G (see §7.2).More directly, one can go back and generalize the basic theorem (7.6.1) to the caseof nonuniform sampling: Suppose that points xi are chosen within the volume V with aprobability density p satisfyingZp dV = 1(7.8.2)The generalized fundamental theorem is that the integral of any function f is estimated, usingN sample points xi , .
. . , xN , bys ZZfhf 2 /p2 i − hf /pi2fI ≡ f dV =±pdV ≈(7.8.3)ppNwhere 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 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 alreadyseen that the idea is to make h = f /p as close to constant as possible. We can be morerigorous by focusing on the numerator inside the square root in equation (7.8.3), which isthe variance per sample point. Both angle brackets are themselves Monte Carlo estimatorsof integrals, so we can write 2 2 Z 2Z2 Z 2Z2fffffS≡≈pdV−=(7.8.4)pdVdV−fdV−p2pp2ppWe now find the optimal p subject to the constraint equation (7.8.2) by the functional variation!Z2ZZ 2fδ(7.8.5)dV −f dV + λ p dV0=δppSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).This section discusses more advanced techniques of Monte Carlo integration. Asexamples 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 wediscuss all fall under the general rubric of reduction of variance (§7.6), but are otherwisequite distinct..
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
