c7-2 (779516)
Текст из файла
2877.2 Transformation Method: Exponential and Normal Deviates7.2 Transformation Method: Exponential andNormal DeviatesThe probability distribution p(x) is of course normalized, so thatZ ∞p(x)dx = 1(7.2.2)−∞Now suppose that we generate a uniform deviate x and then take some prescribedfunction of it, y(x). The probability distribution of y, denoted p(y)dy, is determinedby the fundamental transformation law of probabilities, which is simply|p(y)dy| = |p(x)dx| dx p(y) = p(x) dyor(7.2.3)(7.2.4)Exponential DeviatesAs an example, suppose that y(x) ≡ − ln(x), and that p(x) is as given byequation (7.2.1) for a uniform deviate.
Then dx (7.2.5)p(y)dy = dy = e−y dydywhich is distributed exponentially. This exponential distribution occurs frequentlyin real problems, usually as the distribution of waiting times between independentPoisson-random events, for example the radioactive decay of nuclei. You can alsoeasily see (from 7.2.4) that the quantity y/λ has the probability distribution λe−λy .So we have#include <math.h>float expdev(long *idum)Returns an exponentially distributed, positive,ran1(idum) as the source of uniform deviates.{float ran1(long *idum);float dum;dodum=ran1(idum);while (dum == 0.0);return -log(dum);}random deviate of unit mean,usingSample 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).In the previous section, we learned how to generate random deviates witha uniform probability distribution, so that the probability of generating a numberbetween x and x + dx, denoted p(x)dx, is given byndx 0 < x < 1(7.2.1)p(x)dx =0otherwise288Chapter 7.Random Numbers1uniformdeviate inF(y) =⌡⌠0 p(y)dyyx0ytransformeddeviate outFigure 7.2.1.Transformation method for generating a random deviate y from a known probabilitydistribution p(y). The indefinite integral of p(y) must be known and invertible.
A uniform deviate x ischosen between 0 and 1. Its corresponding y on the definite-integral curve is the desired deviate.Let’s see what is involved in using the above transformation method to generatesome arbitrary desired distribution of y’s, say one with p(y) = f(y) for somepositive function f whose integral is 1. (See Figure 7.2.1.) According to (7.2.4),we need to solve the differential equationdx= f(y)dy(7.2.6)But the solution of this is just x = F (y), where F (y) is the indefinite integral off(y). The desired transformation which takes a uniform deviate into one distributedas f(y) is thereforey(x) = F −1 (x)(7.2.7)where F −1 is the inverse function to F .
Whether (7.2.7) is feasible to implementdepends on whether the inverse function of the integral of f(y) is itself feasible tocompute, either analytically or numerically. Sometimes it is, and sometimes it isn’t.Incidentally, (7.2.7) has an immediate geometric interpretation: Since F (y) isthe area under the probability curve to the left of y, (7.2.7) is just the prescription:choose a uniform random x, then find the value y that has that fraction x ofprobability area to its left, and return the value y.Normal (Gaussian) DeviatesTransformation methods generalize to more than one dimension.
If x1 , x2 ,. . . are random deviates with a joint probability distribution p(x1 , x2 , . . .)dx1 dx2 . . . , and if y1 , y2 , . . . are each functions of all the x’s (same number ofy’s as x’s), then the joint probability distribution of the y’s is ∂(x1 , x2 , . . .) dy1 dy2 . . .(7.2.8)p(y1 , y2 , . . .)dy1 dy2 . .
. = p(x1 , x2 , . . .) ∂(y1 , y2 , . . .) where |∂( )/∂( )| is the Jacobian determinant of the x’s with respect to the y’s(or reciprocal of the Jacobian determinant of the y’s with respect to the x’s).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).p(y)7.2 Transformation Method: Exponential and Normal Deviates289An important example of the use of (7.2.8) is the Box-Muller method forgenerating random deviates with a normal (Gaussian) distribution,21p(y)dy = √ e−y /2 dy2π(7.2.9)Equivalently we can write1x1 = exp − (y12 + y22 )2y21arctanx2 =2πy1Now the Jacobian determinant can readily be calculated (try it!): ∂x∂x11∂(x1 , x2) ∂y1 ∂y2 1 −y12 /21 −y22 /2√ e= ∂x2 ∂x2 = − √ e∂(y1 , y2 )2π2π∂y1(7.2.11)(7.2.12)∂y2Since this is the product of a function of y2 alone and a function of y1 alone, we seethat each y is independently distributed according to the normal distribution (7.2.9).One further trick is useful in applying (7.2.10).
Suppose that, instead of pickinguniform deviates x1 and x2 in the unit square, we instead pick v1 and v2 as theordinate and abscissa of a random point inside the unit circle around the origin. Thenthe sum of their squares, R2 ≡ v12 +v22 is a uniform deviate, which can be used for x1 ,while the angle that (v1 , v2 ) defines with respect to the v1 axis can serve as the randomangle 2πx2 . What’s√ the advantage?√ It’s that the cosine and sine in (7.2.10) can nowbe written as v1 / R2 and v2 / R2 , obviating the trigonometric function calls!We thus have#include <math.h>float gasdev(long *idum)Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum)as the source of uniform deviates.{float ran1(long *idum);static int iset=0;static float gset;float fac,rsq,v1,v2;if (*idum < 0) iset=0;if (iset == 0) {do {v1=2.0*ran1(idum)-1.0;v2=2.0*ran1(idum)-1.0;rsq=v1*v1+v2*v2;Reinitialize.We don’t have an extra deviate handy, sopick two uniform numbers in the square extending from -1 to +1 in each direction,see if they are in the unit circle,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).Consider the transformation between two uniform deviates on (0,1), x1 , x2 , andtwo quantities y1 , y2 ,py1 = −2 ln x1 cos 2πx2(7.2.10)py2 = −2 ln x1 sin 2πx2290Chapter 7.Random Numbers}See Devroye [1] and Bratley [2] for many additional algorithms.CITED REFERENCES AND FURTHER READING:Devroye, L. 1986, Non-Uniform Random Variate Generation (New York: Springer-Verlag), §9.1.[1]Bratley, P., Fox, B.L., and Schrage, E.L.
1983, A Guide to Simulation (New York: SpringerVerlag). [2]Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming(Reading, MA: Addison-Wesley), pp. 116ff.7.3 Rejection Method: Gamma, Poisson,Binomial DeviatesThe rejection method is a powerful, general technique for generating randomdeviates whose distribution function p(x)dx (probability of a value occurring betweenx and x + dx) is known and computable. The rejection method does not requirethat the cumulative distribution function [indefinite integral of p(x)] be readilycomputable, much less the inverse of that function — which was required for thetransformation method in the previous section.The rejection method is based on a simple geometrical argument:Draw a graph of the probability distribution p(x) that you wish to generate, sothat the area under the curve in any range of x corresponds to the desired probabilityof generating an x in that range.
If we had some way of choosing a random point intwo dimensions, with uniform probability in the area under your curve, then the xvalue of that random point would have the desired distribution.Now, on the same graph, draw any other curve f(x) which has finite (notinfinite) area and lies everywhere above your original probability distribution. (Thisis always possible, because your original curve encloses only unit area, by definitionof probability.) We will call this f(x) the comparison function. Imagine nowthat you have some way of choosing a random point in two dimensions that isuniform in the area under the comparison function.
Whenever that point lies outsidethe area under the original probability distribution, we will reject it and chooseanother random point. Whenever it lies inside the area under the original probabilitydistribution, we will accept it. It should be obvious that the accepted points areuniform in the accepted area, so that their x values have the desired distribution. ItSample 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).} while (rsq >= 1.0 || rsq == 0.0);and if they are not, try again.fac=sqrt(-2.0*log(rsq)/rsq);Now make the Box-Muller transformation to get two normal deviates. Return one andsave the other for next time.gset=v1*fac;iset=1;Set flag.return v2*fac;} else {We have an extra deviate handy,iset=0;so unset the flag,return gset;and return it.}.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
