c7-8 (779522), страница 5
Текст из файла (страница 5)
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).float *regn_temp;unsigned long n,npre,nptl,nptr;int j,jb;float avel,varl;float fracl,fval;float rgl,rgm,rgr,s,sigl,siglb,sigr,sigrb;float sum,sumb,summ,summ2;float *fmaxl,*fmaxr,*fminl,*fminr;float *pt,*rmid;7.8 Adaptive and Recursive Monte Carlo Methods327}free_vector(pt,1,ndim);}The miser routine calls a short function ranpt to get a random point within a specifiedd-dimensional region.
The following version of ranpt makes consecutive calls to a uniformrandom number generator and does the obvious scaling. One can easily modify ranpt togenerate its points via the quasi-random routine sobseq (§7.7). We find that miser withsobseq can be considerably more accurate than miser with uniform random deviates. Sincethe use of RSS and the use of quasi-random numbers are completely separable, however, wehave not made the code given here dependent on sobseq.
A similar remark might be maderegarding importance sampling, which could in principle be combined with RSS. (One couldin principle combine vegas and miser, although the programming would be intricate.)extern long idum;void ranpt(float pt[], float regn[], int n)Returns a uniformly random point pt in an n-dimensional rectangular region. Used by miser;calls ran1 for uniform deviates. Your main program should initialize the global variable idumto a negative seed integer.{float ran1(long *idum);int j;for (j=1;j<=n;j++)pt[j]=regn[j]+(regn[n+j]-regn[j])*ran1(&idum);}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).}}free_vector(fminr,1,ndim);free_vector(fminl,1,ndim);free_vector(fmaxr,1,ndim);free_vector(fmaxl,1,ndim);if (!jb) jb=1+(ndim*iran)/175000;MNPT may be too small.rgl=regn[jb];Apportion the remaining points betweenrgm=rmid[jb];left and right.rgr=regn[ndim+jb];fracl=fabs((rgm-rgl)/(rgr-rgl));nptl=(unsigned long)(MNPT+(npts-npre-2*MNPT)*fracl*siglb/(fracl*siglb+(1.0-fracl)*sigrb));Equation (7.8.23).nptr=npts-npre-nptl;regn_temp=vector(1,2*ndim);Now allocate and integrate the two subfor (j=1;j<=ndim;j++) {regions.regn_temp[j]=regn[j];regn_temp[ndim+j]=regn[ndim+j];}regn_temp[ndim+jb]=rmid[jb];miser(func,regn_temp,ndim,nptl,dith,&avel,&varl);regn_temp[jb]=rmid[jb];Dispatch recursive call; will return backregn_temp[ndim+jb]=regn[ndim+jb];here eventually.miser(func,regn_temp,ndim,nptr,dith,ave,var);free_vector(regn_temp,1,2*ndim);*ave=fracl*avel+(1-fracl)*(*ave);*var=fracl*fracl*varl+(1-fracl)*(1-fracl)*(*var);Combine left and right regions by equation (7.8.11) (1st line).free_vector(rmid,1,ndim);Random NumbersChapter 7.328CITED REFERENCES AND FURTHER READING:Hammersley, J.M.
and Handscomb, D.C. 1964, Monte Carlo Methods (London: Methuen).Kalos, M.H. and Whitlock, P.A. 1986, Monte Carlo Methods (New York: Wiley).Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer-Verlag).Press, W.H., and Farrar, G.R. 1990, Computers in Physics, vol. 4, pp. 190–195. [3]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).Lepage, G.P. 1980, “VEGAS: An Adaptive Multidimensional Integration Program,” PublicationCLNS-80/447, Cornell University. [2]Lepage, G.P. 1978, Journal of Computational Physics, vol.
27, pp. 192–203. [1].















