Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 74
Текст из файла (страница 74)
(Compare §7.2.)In general these properties will not hold. Question: What then? Answer: Pull outof the integrand the “best” factor that can be integrated and inverted. The criterionfor “best” is to try to reduce the remaining integrand to a function that is as closeas possible to constant.The limiting case is instructive: If you manage to make the integrand f exactlyconstant, and if the region V , of known volume, exactly encloses the desired regionW , then the average of f that you compute will be exactly its constant value, and theerror estimate in equation (7.6.1) will exactly vanish. You will, in fact, have donethe integral exactly, and the Monte Carlo numerical evaluations are superfluous.
So,backing off from the extreme limiting case, to the extent that you are able to make fapproximately constant by change of variable, and to the extent that you can sample aregion only slightly larger than W , you will increase the accuracy of the Monte Carlointegral. This technique is generically called reduction of variance in the literature.The fundamental disadvantage of simple Monte Carlo integration is that itsaccuracy increases only as the square root of N , the number of sampled points. Ifyour accuracy requirements are modest, or if your computer budget is large, thenthe technique is highly recommended as one of great generality. In the next twosections we will see that there are techniques available for “breaking the square rootof N barrier” and achieving, at least in some cases, higher accuracy with fewerfunction evaluations.7.7 Quasi- (that is, Sub-) Random Sequences309CITED 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 SequencesWe have just seen that choosing N points uniformly randomly in an ndimensionalspace leads to an error term in Monte Carlo integration that decreases√as 1/ N . In essence, each new point sampled adds linearly to an accumulated sumthat will become the function average, and also linearly to an accumulated sum ofsquares that will become the variance (equation 7.6.2). The estimated error comesfrom the square root of this variance, hence the power N −1/2 .Just because this square root convergence is familiar does not, however, meanthat it is inevitable.
A simple counterexample is to choose sample points that lieon a Cartesian grid, and to sample each grid point exactly once (in whatever order).The Monte Carlo method thus becomes a deterministic quadrature scheme — albeita simple one — whose fractional error decreases at least as fast as N −1 (even fasterif the function goes to zero smoothly at the boundaries of the sampled region, oris periodic in the region).The trouble with a grid is that one has to decide in advance how fine it shouldbe.
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 samplepoints “at random,” yet spread out in some self-avoiding way, avoiding the chanceclustering that occurs with uniformly random points.A similar question arises for tasks other than Monte Carlo integration. We mightwant to search an n-dimensional space for a point where some (locally computable)condition holds. Of course, for the task to be computationally meaningful, therehad better be continuity, so that the desired condition will hold in some finite ndimensional neighborhood.
We may not know a priori how large that neighborhoodis, however. We want to “sample until” the desired point is found, moving smoothlyto 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 filln-space more uniformly than uncorrelated random points are called quasi-randomsequences.
That term is somewhat of a misnomer, since there is nothing “random”about quasi-random sequences: They are cleverly crafted to be, in fact, sub-random.The sample points in a quasi-random sequence are, in a precise sense, “maximallyavoiding” of each other.A conceptually simple example is Halton’s sequence [1]. In one dimension, thejth number Hj in the sequence is obtained by the following steps: (i) Write j as anumber in base b, where b is some prime. (For example j = 17 in base b = 3 is122.) (ii) Reverse the digits and put a radix point (i.e., a decimal point base b) in310Chapter 7..
. . . . . . .. . . ... ... . . .. . . ....8.. . . .. . . .. . . . . . . ... .. . .. .. .. ..6.. . . .. . . . . . . . .... . ..4 .. . .. . .. . .. . . .. .... . . . . . . .. . .. . ..2. .. .... .. ........ . . . .010.2.4.6points 1 to 128.8...... .............................. .................. . .. ... . . ............... .......8 ... ... .. ...
. . . . ..... ......................................................... . . . . . . ....6.............................................................. . . ..... . ...4 ....... .. .... ........ ........................ .......... ....... .......................2 .. . ... ... . . .. .. . .. . . .... ......... .................... ....... .............. . .. .0 . .
...2.4.6.8points 513 to 1024...... ..................... .......... ............. . . .... ...... ......... ........8... .. .. . . ........... . . ......... ...... .. .... ........ .... .....6 ...... .. .... ... .. ......... . .... ... ........... .... . ........ . ... ...... ...........4 ... .......... .............. ... ..
.... .. .. . .. . .. . . ..2 . .. . ... . . . . . .. .... . . .......... .... . .... ....... .... ...0 .... .. . . . . .. . .......10110Random Numbers1.2.4.6.8points 129 to 5121.. . ..... . . . ......................................................................................................... .. . . . .. .
. . ... ...8 ............ ...................... ....... .... ............... ................................ ................. ........... ........ ........... .........6 ................................... ........... ...................................................................................... . . . . ..................................4 ... ... ..
............. ......... ................ ................ ...................................................................................2 . .. . ............................................................................................................0 ...
...... .... . .. ....... ... ..10.2.4.6.8points 1 to 10241Figure 7.7.1.First 1024 points of a two-dimensional Sobol’ sequence. The sequence is generatednumber-theoretically, rather than randomly, so successive points at any stage “know” how to fill in thegaps in the previously generated distribution.front of the sequence. (In the example, we get 0.221 base 3.) The result is Hj . Toget a sequence of n-tuples in n-space, you make each component a Halton sequencewith a different prime base b.
Typically, the first n primes are used.It is not hard to see how Halton’s sequence works: Every time the number ofdigits in j increases by one place, j’s digit-reversed fraction becomes a factor ofb finer-meshed. Thus the process is one of filling in all the points on a sequenceof finer and finer Cartesian grids — and in a kind of maximally spread-out orderon each grid (since, e.g., the most rapidly changing digit in j controls the mostsignificant digit of the fraction).Other ways of generating quasi-random sequences have been suggested byFaure, Sobol’, Niederreiter, and others.
Bratley and Fox [2] provide a good reviewand references, and discuss a particularly efficient variant of the Sobol’ [3] sequencesuggested by Antonov and Saleev [4]. It is this Antonov-Saleev variant whoseimplementation we now discuss.3117.7 Quasi- (that is, Sub-) Random SequencesDegreePrimitive Polynomials Modulo 2*10 (i.e., x + 1)21 (i.e., x2 + x + 1)31, 2 (i.e., x3 + x + 1 and x3 + x2 + 1)41, 4 (i.e., x4 + x + 1 and x4 + x3 + 1)52, 4, 7, 11, 13, 1461, 13, 16, 19, 22, 2571, 4, 7, 8, 14, 19, 21, 28, 31, 32, 37, 41, 42, 50, 55, 56, 59, 62814, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84, 97, 103, 115, 12298, 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, 253104, 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 thehigh-order bit and the unit bit).The Sobol’ sequencegenerates numbers between zero and one directly as binary fractionsof length w bits, from a set of w special binary fractions, Vi , i = 1, 2, .
. . , w, called directionnumbers. In Sobol’s original method, the jth number Xj is generated by XORing (bitwiseexclusive or) together the set of Vi ’s satisfying the criterion on i, “the ith bit of j is nonzero.”As j increments, in other words, different ones of the Vi ’s flash in and out of Xj on differenttime scales. V1 alternates between being present and absent most quickly, while Vk goes frompresent to absent (or vice versa) only every 2k−1 steps.Antonov and Saleev’s contribution was to show that instead of using the bits of theinteger j to 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.)Now G(j) and G(j + 1) differ in exactly one bit position, namely in the position of therightmost zero bit in the binary representation of j (adding a leading zero to j if necessary).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
