Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 75
Текст из файла (страница 75)
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)312Chapter 7.Random NumbersInitializing Values Used in sobseqDegreePolynomialStarting 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.7.7 Quasi- (that is, Sub-) Random Sequences313in=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.}}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 by$1/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, giving& & &dx 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’314Chapter 7.Random Numbersfractional accuracy of integral.1∝ 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”). 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,01r < 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.7.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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
