c7-7 (779521), страница 2
Текст из файла (страница 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). Aconsequence is that the j + 1st Sobol’-Antonov-Saleev number can be obtained from the jthby XORing it with a single Vi , namely with i the position of the rightmost zero bit in j.
Thismakes the calculation of the sequence very efficient, as we shall see.Figure 7.7.1 plots the first 1024 points generated by a two-dimensional Sobol’ sequence.One sees that successive points do “know” about the gaps left previously, and keep fillingthem in, hierarchically.We have deferred to this point a discussion of how the direction numbers Vi are generated.Some nontrivial mathematics is involved in that, so we will content ourself with a cookbooksummary only: Each different Sobol’ sequence (or component of an n-dimensional sequence)is based on a different primitive polynomial over the integers modulo 2, that is, a polynomialwhose coefficients are either 0 or 1, and which generates a maximal length shift registersequence. (Primitive polynomials modulo 2 were used in §7.4, and are further discussed in§20.3.) Suppose P is such a polynomial, of degree q,P = xq + a1 xq−1 + a2 xq−2 + · · · + aq−1 + 1(7.7.1)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).1312Chapter 7.Random NumbersInitializing Values Used in sobseqPolynomialStarting Values101(3)(5)(15) .
. .2111(7)(11) . . .31137(5) . . .32133(15) . . .4111313 . . .441159 ...Parenthesized values are not freely specifiable, but are forced by the required recurrencefor this degree.Define a sequence of integers Mi by the q-term recurrence relation,Mi = 2a1 Mi−1 ⊕ 22 a2 Mi−2 ⊕ · · · ⊕ 2q−1 Mi−q+1 aq−1 ⊕ (2q Mi−q ⊕ Mi−q ) (7.7.2)Here bitwise XOR is denoted by ⊕. The starting values for this recurrence are that M1, . . .
, Mqcan be arbitrary odd integers less than 2, . . . , 2q , respectively. Then, the direction numbersVi are given byVi = Mi /2ii = 1, . . . , w(7.7.3)The accompanying table lists all primitive polynomials modulo 2 with degree q ≤ 10.Since the coefficients are either 0 or 1, and since the coefficients of xq and of 1 are predictably1, it is convenient to denote a polynomial by its middle coefficients taken as the bits of a binarynumber (higher powers of x being more significant bits). The table uses this convention.Turn now to the implementation of the Sobol’ sequence. Successive calls to the functionsobseq (after a preliminary initializing call) return successive points in an n-dimensionalSobol’ sequence based on the first n primitive polynomials in the table.
As given, theroutine is initialized for maximum n of 6 dimensions, and for a word length w of 30 bits.These parameters can be altered by changing MAXBIT (≡ w) and MAXDIM, and by addingmore initializing data to the arrays ip (the primitive polynomials from the table), mdeg (theirdegrees), and iv (the starting values for the recurrence, equation 7.7.2). A second table,above, elucidates the initializing data in the routine.#include "nrutil.h"#define MAXBIT 30#define MAXDIM 6void sobseq(int *n, float x[])When n is negative, internally initializes a set of MAXBIT direction numbers for each of MAXDIMdifferent Sobol’ sequences.
When n is positive (but ≤MAXDIM), returns as the vector x[1..n]the next values from n of these sequences. (n must not be changed between initializations.){int j,k,l;unsigned long i,im,ipp;static float fac;static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};static unsigned long iv[MAXDIM*MAXBIT+1]={0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};if (*n < 0) {for (k=1;k<=MAXDIM;k++) ix[k]=0;Initialize, don’t return a vector.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).Degree7.7 Quasi- (that is, Sub-) Random Sequences313}How good is a Sobol’ sequence, anyway? For Monte Carlo integration of a smoothfunction in n dimensions, the answer is that the fractional error will decrease with N , thenumber of samples, as (ln N )n /N , i.e., almost as fast as 1/N .
As an example, let us integratea function that is nonzero inside a torus (doughnut) in three-dimensional space. If the majorradius of the torus is R0 , the minor radial coordinate r is defined by1/2(7.7.4)r = [(x2 + y2 )1/2 − R0 ]2 + z 2Let us try the functionf (x, y, z) =1 + cos0πr2a2r < r0(7.7.5)r ≥ r0which can be integrated analytically in cylindrical coordinates, givingZ Z Zdx dy dz f (x, y, z) = 2π2 a2 R0(7.7.6)With parameters R0 = 0.6, r0 = 0.3, we did 100 successive Monte Carlo integrations ofequation (7.7.4), sampling uniformly in the region −1 < x, y, z < 1, for the two cases ofuncorrelated random points and the Sobol’ sequence generated by the routine sobseq. Figure7.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 single integration, the error of course wandersfrom 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 showsthe familiar N −1/2 asymptotics. The thin, solid gray curve shows the result for the Sobol’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).in=0;if (iv[1] != 1) return;fac=1.0/(1L << MAXBIT);for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];To allow both 1D and 2D addressing.for (k=1;k<=MAXDIM;k++) {for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j);Stored values only require normalization.for (j=mdeg[k]+1;j<=MAXBIT;j++) {Use the recurrence to get other valipp=ip[k];ues.i=iu[j-mdeg[k]][k];i ^= (i >> mdeg[k]);for (l=mdeg[k]-1;l>=1;l--) {if (ipp & 1) i ^= iu[j-l][k];ipp >>= 1;}iu[j][k]=i;}}} else {Calculate the next vector in the seim=in++;quence.for (j=1;j<=MAXBIT;j++) {Find the rightmost zero bit.if (!(im & 1)) break;im >>= 1;}if (j > MAXBIT) nrerror("MAXBIT too small in sobseq");im=(j-1)*MAXDIM;for (k=1;k<=IMIN(*n,MAXDIM);k++) {XOR the appropriate direction numix[k] ^= iv[im+k];ber into each component of thex[k]=ix[k]*fac;vector and convert to a floating}number.}314Chapter 7.Random Numbers∝ N −1/2.01.001100∝ N −2/3pseudo-random, hard boundarypseudo-random, soft boundaryquasi-random, hard boundaryquasi-random, soft boundary1000∝ N −11000010 5number of points NFigure 7.7.2.
Fractional accuracy of Monte Carlo integrations as a function of number of points 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. Quasirandom sampling does better when the integrand is smooth (“soft boundary”) than when it has stepdiscontinuities (“hard boundary”).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
