7.8 Adaptive and Recursive Monte Carlo Methods 325 find empirically that it is somewhat more robust to use the square of the difference of maximum and minimum sampled function values,instead of the genuine second moment of the samples. This estimator is of course increasingly biased with increasing sample size;however,equation (7.8.23)uses it only to compare two subvolumes(a and b)having approximately equal numbers of samples.The "max minus min"estimator proves its worth when the preliminary sampling yields only a single point,or small number of points,in active regions of the integrand.In many realistic cases,these are indicators of nearby regions of even 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 ofa"dithering parameter," dith,whose nonzero value causes subvolumes to be divided not exactly down the middle,but rather into fractions 0.5+dith,with the sign of the randomly chosen by a built-in random number routine.Normally dith can be set to zero.However,there is a large advantage in 8 taking dith to be nonzero if some special symmetry of the integrand puts the active region exactly at the midpoint of the region,or at the center of some power-of-two submultiple of the region.One wants to avoid the extreme case of the active region being evenly divided nted for into 2d abutting corners of a d-dimensional space.A typical nonzero value of dith,on those occasions when it is useful,might be 0.1.Of course,when the dithering parameter is nonzero,we must take the differing sizes of the subvolumes into account;the code does this through the variable fracl. One final feature in the code deserves mention.The RSS algorithm uses a single set of sample points to evaluate equation (7.8.23)in all d directions.At bottom levels of the recursion,the number of sample points can be quite small.Although rare,it can happen that in one direction all the samples are in one half of the volume;in that case,that direction is ignored as a candidate for bifurcation.Even more rare is the possibility that all of the samples are in one half of the volume in al/directions.In this case,a random direction is 号&s3d的 Press. ART chosen.If this happens too often in your application,then you should increase MNPT(see line if (!ib)...in the code). Note that miser,as given,returns as ave an estimate of the average function value 9 Programs (f)),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 SCIENTIFIC (7.8.8),since the volume V of the rectangular region is known. #include <stdlib.h> #include <math.h> COMPUTING #include "nrutil.h" 共daf3 ne peac0.1 19881992 #define MNPT 15 (ISBN #define MNBS 60 #define TINY 1.0e-30 色 #define BIG 1.0e30 Here PFAC is the fraction of remaining function evaluations used at each stage to explore the variance of func.At least MNPT function evaluations are performed in any terminal subregion; Numerical Recipes 43108 a subregion is further bisected only if at least MNBS function evaluations are available.We take MNBS =4*MNPT. (outside static long iran=0; North Software. 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 volume specified by regn[1..2*ndim],a vector consisting of ndim "lower-left"coordinates of the region followed by ndim "upper-right"coordinates.The function is sampled a total of npts times,at locations determined by the method of recursive stratified sampling.The mean value of 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 normally be set to zero,but can be set to (e.g.)0.1 if func's active region falls on the boundary of a power-of-two subdivision of region void ranpt(float pt[],float regn[],int n); float *regn_temp;7.8 Adaptive and Recursive Monte Carlo Methods 325 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). find empirically that it is somewhat more robust to use the square of the difference of maximum and minimum sampled function values, instead of the genuine second moment of the samples. This estimator is of course increasingly biased with increasing sample size; however, equation (7.8.23) uses it only to compare two subvolumes (a and b) having approximately equal numbers of samples. The “max minus min” estimator proves its worth when the preliminary sampling yields only a single point, or small number of points, in active regions of the integrand. In many realistic cases, these are indicators of nearby regions of even 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, but rather into fractions 0.5±dith, with the sign of the ± randomly chosen by a built-in random number routine. Normally dith can be set to zero. However, there is a large advantage in taking dith to be nonzero if some special symmetry of the integrand puts the active region exactly at the midpoint of the region, or at the center of some power-of-two submultiple of the region. One wants to avoid the extreme case of the active region being evenly divided into 2d abutting corners of a d-dimensional space. A typical nonzero value of dith, on those occasions when it is useful, might be 0.1. Of course, when the dithering parameter is nonzero, we must take the differing sizes of the subvolumes into account; the code does this through the variable fracl. One final feature in the code deserves mention. The RSS algorithm uses a single set of sample points to evaluate equation (7.8.23) in all d directions. At bottom levels of the recursion, the number of sample points can be quite small. Although rare, it can happen that in one direction all the samples are in one half of the volume; in that case, that direction is ignored as a candidate for bifurcation. Even more rare is the possibility that all of the samples are in one half of the volume in all directions. In this case, a random direction is chosen. If this happens too often in your application, then you should increase MNPT (see line if (!jb)... in the code). Note that miser, as given, returns as ave an estimate of the average function value f, 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. #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.0e30 Here PFAC is the fraction of remaining function evaluations used at each stage to explore the variance 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 take MNBS = 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 volume specified by regn[1..2*ndim], a vector consisting of ndim “lower-left” coordinates of the region followed by ndim “upper-right” coordinates. The function is sampled a total of npts times, at locations determined by the method of recursive stratified sampling. The mean value of 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 normally be set to zero, but can be set to (e.g.) 0.1 if func’s active region falls on the boundary of a power-of-two subdivision of region. { void ranpt(float pt[], float regn[], int n); float *regn_temp;