Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 78
Текст из файла (страница 78)
Enter here on a cold start.mds=ndo=1;Change to mds=0 to disable stratified sampling,for (j=1;j<=ndim;j++) xi[j][1]=1.0;i.e., use importance sampling only.}if (init <= 1) si=swgt=schi=0.0;Enter here to inherit the grid from a previous call, but not its answers.if (init <= 2) {Enter here to inherit the previous grid and itsnd=NDMX;answers.ng=1;if (mds) {Set up for stratification.ng=(int)pow(ncall/2.0+0.25,1.0/ndim);mds=1;if ((2*ng-NDMX) >= 0) {mds = -1;npg=ng/NDMX+1;nd=ng/npg;ng=npg*nd;}}for (k=1,i=1;i<=ndim;i++) k *= ng;npg=IMAX(ncall/k,2);calls=(float)npg * (float)k;dxg=1.0/ng;for (dv2g=1,i=1;i<=ndim;i++) dv2g *= dxg;dv2g=SQR(calls*dv2g)/npg/npg/(npg-1.0);xnd=nd;dxg *= xnd;xjac=1.0/calls;for (j=1;j<=ndim;j++) {dx[j]=regn[j+ndim]-regn[j];xjac *= dx[j];}if (nd != ndo) {Do binning if necessary.for (i=1;i<=IMAX(nd,ndo);i++) r[i]=1.0;for (j=1;j<=ndim;j++) rebin(ndo/xnd,nd,r,xin,xi[j]);ndo=nd;}if (nprn >= 0) {printf("%s: ndim= %3d ncall= %8.0f\n"," Input parameters for vegas",ndim,calls);printf("%28s it=%5d itmx=%5d\n"," ",it,itmx);printf("%28s nprn=%3d ALPH=%5.2f\n"," ",nprn,ALPH);printf("%28s mds=%3d nd=%4d\n"," ",mds,nd);for (j=1;j<=ndim;j++) {printf("%30s xl[%2d]= %11.4g xu[%2d]= %11.4g\n"," ",j,regn[j],j,regn[j+ndim]);}}}for (it=1;it<=itmx;it++) {Main iteration loop.
Can enter here (init ≥ 3) to do an additional itmx iterations withall other parameters unchanged.ti=tsi=0.0;for (j=1;j<=ndim;j++) {kg[j]=1;for (i=1;i<=nd;i++) d[i][j]=di[i][j]=0.0;322Chapter 7.Random Numbers}for (;;) {fb=f2b=0.0;for (k=1;k<=npg;k++) {wgt=xjac;for (j=1;j<=ndim;j++) {xn=(kg[j]-ran2(&idum))*dxg+1.0;ia[j]=IMAX(IMIN((int)(xn),NDMX),1);if (ia[j] > 1) {xo=xi[j][ia[j]]-xi[j][ia[j]-1];rc=xi[j][ia[j]-1]+(xn-ia[j])*xo;} else {xo=xi[j][ia[j]];rc=(xn-ia[j])*xo;}x[j]=regn[j]+rc*dx[j];wgt *= xo*xnd;}f=wgt*(*fxn)(x,wgt);f2=f*f;fb += f;f2b += f2;for (j=1;j<=ndim;j++) {di[ia[j]][j] += f;if (mds >= 0) d[ia[j]][j] += f2;}}f2b=sqrt(f2b*npg);f2b=(f2b-fb)*(f2b+fb);if (f2b <= 0.0) f2b=TINY;ti += fb;tsi += f2b;if (mds < 0) {Use stratified sampling.for (j=1;j<=ndim;j++) d[ia[j]][j] += f2b;}for (k=ndim;k>=1;k--) {kg[k] %= ng;if (++kg[k] != 1) break;}if (k < 1) break;}tsi *= dv2g;Compute final results for this iteration.wgt=1.0/tsi;si += wgt*ti;schi += wgt*ti*ti;swgt += wgt;*tgral=si/swgt;*chi2a=(schi-si*(*tgral))/(it-0.9999);if (*chi2a < 0.0) *chi2a = 0.0;*sd=sqrt(1.0/swgt);tsi=sqrt(tsi);if (nprn >= 0) {printf("%s %3d : integral = %14.7g +/- %9.2g\n"," iteration no.",it,ti,tsi);printf("%s integral =%14.7g+/-%9.2g chi**2/IT n = %9.2g\n"," all iterations: ",*tgral,*sd,*chi2a);if (nprn) {for (j=1;j<=ndim;j++) {printf(" DATA FOR axis %2d\n",j);printf("%6s%13s%11s%13s%11s%13s\n","X","delta i","X","delta i","X","delta i");for (i=1+nprn/2;i<=nd;i += nprn+2) {printf("%8.5f%12.4g%12.5f%12.4g%12.5f%12.4g\n",xi[j][i],di[i][j],xi[j][i+1],7.8 Adaptive and Recursive Monte Carlo Methods323di[i+1][j],xi[j][i+2],di[i+2][j]);}}}}for (j=1;j<=ndim;j++) {Refine the grid.
Consult references to understandxo=d[1][j];the subtlety of this procedure. The refinexn=d[2][j];ment is damped, to avoid rapid, destabilizd[1][j]=(xo+xn)/2.0;ing changes, and also compressed in rangedt[j]=d[1][j];by the exponent ALPH.for (i=2;i<nd;i++) {rc=xo+xn;xo=xn;xn=d[i+1][j];d[i][j] = (rc+xn)/3.0;dt[j] += d[i][j];}d[nd][j]=(xo+xn)/2.0;dt[j] += d[nd][j];}for (j=1;j<=ndim;j++) {rc=0.0;for (i=1;i<=nd;i++) {if (d[i][j] < TINY) d[i][j]=TINY;r[i]=pow((1.0-d[i][j]/dt[j])/(log(dt[j])-log(d[i][j])),ALPH);rc += r[i];}rebin(rc/xnd,nd,r,xin,xi[j]);}}}void rebin(float rc, int nd, float r[], float xin[], float xi[])Utility routine used by vegas, to rebin a vector of densities xi into new bins defined by avector r.{int i,k=0;float dr=0.0,xn=0.0,xo=0.0;for (i=1;i<nd;i++) {while (rc > dr)dr += r[++k];if (k > 1) xo=xi[k-1];xn=xi[k];dr -= rc;xin[i]=xn-(xn-xo)*dr/r[k];}for (i=1;i<nd;i++) xi[i]=xin[i];xi[nd]=1.0;}Recursive Stratified SamplingThe problem with stratified sampling, we have seen, is that it may not avoid the K dexplosion inherent in the obvious, Cartesian, tesselation of a d-dimensional volume.
Atechnique called recursive stratified sampling [3] attempts to do this by successive bisectionsof a volume, not along all d dimensions, but rather along only one dimension at a time.324Chapter 7.Random NumbersThe starting points are equations (7.8.10) and (7.8.13), applied to bisections of successivelysmaller subregions.Suppose that we have a quota of N evaluations of the function f , and want to evaluatef in the rectangular parallelepiped region R = (xa , xb ). (We denote such a region by thetwo coordinate vectors of its diagonally opposite corners.) First, we allocate a fraction p ofN towards exploring the variance of f in R: We sample pN function values uniformly inR and accumulate the sums that will give the d different pairs of variances corresponding tothe d different coordinate directions along which R can be bisected.
In other words, in pNsamples, we estimate Var (f ) in each of the regions resulting from a possible bisection of R,Rai ≡(xa , xb −1ei · (xb − xa )ei )2(7.8.21)1Rbi ≡(xa + ei · (xb − xa )ei , xb )2Here ei is the unit vector in the ith coordinate direction, i = 1, 2, . . . , d.Second, we inspect the variances to find the most favorable dimension i to bisect. Byequation (7.8.15), we could, for example, choose that i for which the sum of the square rootsof the variance estimators in regions Rai and Rbi is minimized. (Actually, as we will explain,we do something slightly different.)Third, we allocate the remaining (1 − p)N function evaluations between the regionsRai and Rbi . If we used equation (7.8.15) to choose i, we should do this allocation accordingto equation (7.8.14).We now have two parallelepipeds each with its own allocation of function evaluationsfor estimating the mean of f .
Our “RSS” algorithm now shows itself to be recursive: Toevaluate the mean in each region, we go back to the sentence beginning “First,...” in theparagraph above equation (7.8.21). (Of course, when the allocation of points to a region fallsbelow some number, we resort to simple Monte Carlo rather than continue with the recursion.)Finally, we combine the means, and also estimated variances of the two subvolumes,using equation (7.8.10) and the first line of equation (7.8.11).This completes the RSS algorithm in its simplest form.
Before we describe someadditional tricks under the general rubric of “implementation details,” we need to returnbriefly to equations (7.8.13)–(7.8.15) and derive the equations that we actually use instead ofthese. The right-hand side of equation (7.8.13) applies the familiar scaling law of equation(7.8.9) twice, once to a and again to b. This would be correct if the estimates f a and f bwere each made by simple Monte Carlo, with uniformly random sample points. However, thetwo estimates of the mean are in fact made recursively.
Thus, there is no reason to expectequation (7.8.9) to hold. Rather, we might substitute for equation (7.8.13) the relation,,1 Vara (f )Varb (f )Var f =+(7.8.22)4Naα(N − Na )αwhere α is an unknown constant ≥ 1 (the case of equality,-corresponding to simple MonteCarlo). In that case, a short calculation shows that Var f is minimized whenNaVara (f )1/(1+α)=1/(1+α)NVara (f )+ Varb (f )1/(1+α)and that its minimum value is1+α- ,Var f ∝ Vara (f )1/(1+α) + Varb (f )1/(1+α)(7.8.23)(7.8.24)Equations (7.8.22)–(7.8.24) reduce to equations (7.8.13)–(7.8.15) when α = 1.
Numericalexperiments to find a self-consistent value for α find that α ≈ 2. That is, when equation(7.8.23) with α = 2 is used recursively to allocate sample opportunities, the observed varianceof the RSS algorithm goes approximately as N −2 , while any other value of α in equation(7.8.23) gives a poorer fall-off. (The sensitivity to α is, however, not very great; it is notknown whether α = 2 is an analytically justifiable result, or only a useful heuristic.)The principal difference between miser’s implementation and the algorithm as describedthus far lies in how the variances on the right-hand side of equation (7.8.23) are estimated.7.8 Adaptive and Recursive Monte Carlo Methods325We find empirically that it is somewhat more robust to use the square of the difference ofmaximum and minimum sampled function values, instead of the genuine second momentof the samples.
This estimator is of course increasingly biased with increasing samplesize; however, equation (7.8.23) uses it only to compare two subvolumes (a and b) havingapproximately equal numbers of samples. The “max minus min” estimator proves its worthwhen the preliminary sampling yields only a single point, or small number of points, in activeregions of the integrand. In many realistic cases, these are indicators of nearby regions ofeven greater importance, and it is useful to let them attract the greater sampling weight that“max minus min” provides.A second modification embodied in the code is the introduction of a “dithering parameter,”dith, whose nonzero value causes subvolumes to be divided not exactly down the middle, butrather into fractions 0.5±dith, with the sign of the ± randomly chosen by a built-in randomnumber routine. Normally dith can be set to zero.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
