c7-8 (779522), страница 4
Текст из файла (страница 4)
. . , 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 hf ia and hf ibwere 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 hf i0 =+(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 hf i0 is minimized whenNaVara (f )1/(1+α)=NVara (f )1/(1+α) + Varb (f )1/(1+α)and that its minimum value isi1+α hVar hf i0 ∝ 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.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 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 evaluatehf i0 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,7.8 Adaptive and Recursive Monte Carlo Methods325#include <stdlib.h>#include <math.h>#include "nrutil.h"#define PFAC 0.1#define MNPT 15#define MNBS 60#define TINY 1.0e-30#define BIG 1.0e30Here PFAC is the fraction of remaining function evaluations used at each stage to explore thevariance of func.
At least MNPT function evaluations are performed in any terminal subregion;a subregion is further bisected only if at least MNBS function evaluations are available. We takeMNBS = 4*MNPT.static long iran=0;void miser(float (*func)(float []), float regn[], int ndim, unsigned long npts,float dith, float *ave, float *var)Monte Carlo samples a user-supplied ndim-dimensional function func in a rectangular volumespecified by regn[1..2*ndim], a vector consisting of ndim “lower-left” coordinates of theregion followed by ndim “upper-right” coordinates.
The function is sampled a total of nptstimes, at locations determined by the method of recursive stratified sampling. The mean valueof the function in the region is returned as ave; an estimate of the statistical uncertainty of ave(square of standard deviation) is returned as var. The input parameter dith should normallybe set to zero, but can be set to (e.g.) 0.1 if func’s active region falls on the boundary of apower-of-two subdivision of region.{void ranpt(float pt[], float regn[], int n);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).We 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. However, there is a large advantage intaking dith to be nonzero if some special symmetry of the integrand puts the active regionexactly at the midpoint of the region, or at the center of some power-of-two submultiple ofthe region.
One wants to avoid the extreme case of the active region being evenly dividedinto 2d abutting corners of a d-dimensional space. A typical nonzero value of dith, onthose occasions when it is useful, might be 0.1. Of course, when the dithering parameteris nonzero, we must take the differing sizes of the subvolumes into account; the code doesthis through the variable fracl.One final feature in the code deserves mention.
The RSS algorithm uses a single setof sample points to evaluate equation (7.8.23) in all d directions. At bottom levels of therecursion, the number of sample points can be quite small. Although rare, it can happen thatin one direction all the samples are in one half of the volume; in that case, that directionis ignored as a candidate for bifurcation. Even more rare is the possibility that all of thesamples are in one half of the volume in all directions. In this case, a random direction ischosen. If this happens too often in your application, then you should increase MNPT (seeline if (!jb).
. . in the code).Note that miser, as given, returns as ave an estimate of the average function valuehhf ii, not the integral of f over the region. The routine vegas, adopting the other convention,returns as tgral the integral.
The two conventions are of course trivially related, by equation(7.8.8), since the volume V of the rectangular region is known.326Chapter 7.Random Numberspt=vector(1,ndim);if (npts < MNBS) {Too few points to bisect; do straightsumm=summ2=0.0;Monte Carlo.for (n=1;n<=npts;n++) {ranpt(pt,regn,ndim);fval=(*func)(pt);summ += fval;summ2 += fval * fval;}*ave=summ/npts;*var=FMAX(TINY,(summ2-summ*summ/npts)/(npts*npts));}else {Do the preliminary (uniform) sampling.rmid=vector(1,ndim);npre=LMAX((unsigned long)(npts*PFAC),MNPT);fmaxl=vector(1,ndim);fmaxr=vector(1,ndim);fminl=vector(1,ndim);fminr=vector(1,ndim);for (j=1;j<=ndim;j++) {Initialize the left and right bounds foriran=(iran*2661+36979) % 175000;each dimension.s=SIGN(dith,(float)(iran-87500));rmid[j]=(0.5+s)*regn[j]+(0.5-s)*regn[ndim+j];fminl[j]=fminr[j]=BIG;fmaxl[j]=fmaxr[j] = -BIG;}for (n=1;n<=npre;n++) {Loop over the points in the sample.ranpt(pt,regn,ndim);fval=(*func)(pt);for (j=1;j<=ndim;j++) {Find the left and right bounds for eachif (pt[j]<=rmid[j]) {dimension.fminl[j]=FMIN(fminl[j],fval);fmaxl[j]=FMAX(fmaxl[j],fval);}else {fminr[j]=FMIN(fminr[j],fval);fmaxr[j]=FMAX(fmaxr[j],fval);}}}sumb=BIG;Choose which dimension jb to bisect.jb=0;siglb=sigrb=1.0;for (j=1;j<=ndim;j++) {if (fmaxl[j] > fminl[j] && fmaxr[j] > fminr[j]) {sigl=FMAX(TINY,pow(fmaxl[j]-fminl[j],2.0/3.0));sigr=FMAX(TINY,pow(fmaxr[j]-fminr[j],2.0/3.0));sum=sigl+sigr;Equation (7.8.24), see text.if (sum<=sumb) {sumb=sum;jb=j;siglb=sigl;sigrb=sigr;}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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
