SCHOOL OF MATHEMATICS PRINCETON, NEW JERSEY ANOTHER VON NEUMANN LETTER Fig. 3. In this 1947 letter to Stan Ulam, von Neumann discusses two methods for generating the nonuniform distributions of random numbers needed in the Monte Carlo method. The second paragraph summarizes the inverse-function approach in which (x') represents the uniform distribution and (g′) the desired nonuniform distribution. The rest of the letter describes an alternative approach based on two uniform and independent distributions: (x′) and (y'). In this latter approach a value x' from the first set is accepted when a value y' from the second set satisfies the condition y' ≤ f(x'), where f(') dε is the density of the desired distribution function. (It should be noted that in von Neumann's example for forming the random pairs & = sin x and ŋ = cos x, he probably meant to say that x is equidistributed between 0 and 360 degrees (rather than “300”). Also, his notation for the tangent function is "tg," so that the second set of equations for έ and ʼn are just half-angle (y = x/2) trigonometric identities.) Put या n عمر = = 2 t with y (which is ) equidistributed between 0° and 180°. 2 1+t 1- t2 7 1+t2 Then the §, will have to be replaced randomly by 7, § This can be done by using random digits 0, } -21 -} 8 Reject this digit lies between 0 and 1, and its distriHence one may pick pairs of numbers t, s both (independently) equidistributed between 0 and 1, and then use t t, s t and at по step } for (1+ t2) s≤ 1 reject form this for (1+ t2)s > I Of course, the first part requires a divider, but the method may still be worth keeping in mind, especially when the ENIAC is available. proposed by Stan, uses the inverse of the desired function f = g. For example, to get the exponentially decreasing distribution of random numbers on the interval (0,∞) needed for path lengths, one applies the inverse function f (x) = − lnx to a uniform distribution of random numbers on the open interval (0, 1). What if it is difficult or computationally expensive to form the inverse function, which is frequently true when the desired function is empirical? The rest of von Neumann's letter describes an alternative technique that will work for such cases. In this approach two uniform and independent distributions (x') and (y') are used. A value x' from the first set is accepted when a value y' from the second set satisfies the condition y' <ƒ(x'), where f()d is the density of the desired distribution function (that is, g(x) = ff(x)dx). This acceptance-rejection technique of von Neumann's can best be illustrated graphically (Fig. 4). If the two numbers x' and y' are selected randomly from the domain and range, respectively, of the function f, then each pair of numbers represents a point in the function's coordinate plane (x', y'). When y' > f(x') the point lies above the curve for f(x), and x1 is rejected; when y' <f(x') the point lies on or below the curve, and x is accepted. Thus, the fraction of accepted points is equal to the fraction of the area below the curve. In fact, the proportion of points selected that fall in a small interval along the x-axis will be proportional to the average height of the curve in that interval, ensuring generation of random numbers that mirror the desired distribution. fter series of "games" have been Aplayed, how does one extract mean ingful information? For each of thousands of neutrons, the variables describing the chain of events are stored, and this collection constitutes a numerical model of the process being studied. The collection of variables is analyzed using sta tistical methods identical to those used to analyze experimental observations of physical processes. One can thus extract information about any variable that was accounted for in the process. For example, the average energy of the neutrons at a particular time is calculated by simply taking the average of all the values generated by the chains at that time. This value has an uncertainty proportional to √V/(N-1), where V is the variance of, in this case, the energy and N is the number of trials, or chains, followed. It is, of course, desirable to reduce statistical uncertainty. Any modification to the stochastic calculational process that generates the same expectation values but smaller variances is called a variance frequently reflect the addition of known. physics to the problem, and they reduce the variance by effectively increasing the number of data points pertinent to the variable of interest. An example is dealing with neutron absorption by weighted sampling. In this technique, each neutron is assigned a unit "weight" at the start of its path. The weight is then decreased, bit by bit at each collision, in proportion to the absorption cross section divided by the total collision cross section. After each collision an outcome other than absorption is selected by random sampling and the path is continued. This technique reduces the variance by replacing the sudden, one-time process of neutron absorption by a gradual elimination of the neutron. Another example of variance reduction is a technique that deals with outcomes that terminate a chain. Say that at each collision one of the alternative outcomes terminates the chain and associated with this outcome is a particular value x, for the variable of interest (an example is x, being a path length long enough for the neutron to escape). Instead of collecting these values only when the chain terminates, one can generate considerably more data about this particular outcome by making an extra calculation at each decision point. In this calculation the know value x, for termination is multiplied by the probability that that outcome will occur. Then random values are selected to continue the chain in the usual manner. By the end of the calculation, the "weighted values" for the terminating outcome have been summed over all decision points. This variance-reduction technique is especially useful if the probablity of the alternative in question is low. For example, shielding calculations typically predict that only one in many thousands of neutrons actually get through the shielding. Instead of accumulating those rare paths, the small probabilities that a neutron will get through the shield on its very next free flight are accumulated after each collision. he Monte Carlo method has proven The fact, "solitaire games" now range from Random-Number R by Tony Warnock andom numbers have applications in many areas: simulation, game-playing, cryptography, statistical sampling, evaluation of multiple integrals, particletransport calculations, and computations in statistical physics, to name a few. Since each application involves slightly different criteria for judging the "worthiness" of the random numbers generated, a variety of generators have been developed, each with its own set of advantages and disadvantages. Depending on the application, three types of number sequences might prove adequate as the "random numbers." From a purist point of view, of course, a series of numbers generated by a truly random process is most desirable. This type of sequence is called a random-number sequence, and one of the key problems is deciding whether or not the generating process is, in fact, random. A more practical sequence is the pseudo-random sequence, a series of numbers generated by a deterministic process that is intended merely to imitate a random sequence but which, of course, does not rigorously obey such things as the laws of large numbers (see page 69). Finally, a quasi-random sequence is a series of numbers that makes no pretense at being random but that has important predefined statistical properties shared with random sequences. Physical Random-Number Generators Games of chance are the classic examples of random processes, and the first inclination would be to use traditional gambling devices as random-number generators. Unfortunately, these devices are rather slow, especially since the typical computer application may require millions of numbers per second. Also, the numbers obtained from such devices are not always truly random: cards may be imperfectly shuffled, dice may not be true, wheels may not be balanced, and so forth. However, in the early 1950s the Rand Corporation constructed a million-digit table of random numbers using an electrical "roulette wheel." (The device had 32 slots, of which 12 were ignored; the others were numbered from 0 to 9 twice.) Classical gambling devices appear random only because of our ignorance of initial conditions; in principle, these devices follow deterministic Newtonian physics. Another possibility for generating truly random numbers is to take advantage of the Heisenberg uncertainty principle and quantum effects, say by counting decays of a radioactive source or by tapping into electrical noise. Both of these methods have been used to generate random numbers for computers, but both suffer the defects of slowness and ill-defined distributions (however, on a different but better order of magnitude than gambling devices). For instance, although each decay in a radioactive source may occur randomly and independently of other decays, it is not necessarily true that successive counts in the detector are independent of each other. The time it takes to reset the counter, for example, might depend on the previous count. Furthermore, the source itself constantly changes in time as the number of remaining radioactive particles decreases exponentially. Also, voltage drifts can introduce bias into the noise of electrical devices. There are, of course, various tricks to overcome some of these disadvantages. One can partially compensate for the counter-reset problem by replacing the string of bits that represents a given count with a new number in which all of the original 1-1 and 0-0 pairs have been discarded and all of the original 0-1 and 1-0 pairs have been changed to 0 and 1, respectively. This trick reduces the bias caused when the probability of a O is different from that of a 1 but does not completely eliminate nonindependence of successive counts. A shortcoming of any physical generator is the lack of reproducibility. Reproducibility is needed for debugging codes that use the random numbers and for making correlated or anti-correlated computations. Of course, if one wants random numbers for a cryptographic one-time pad, reproducibility is the last attribute desired, and time can be traded for security. A radioactive source used with the bias-removal technique described above is probably sufficient. Arithmetical Pseudo-Random Generators The most common method of generating pseudo-random numbers on the computer uses a recursive technique called the linear-congruential, or Lehmer, generator. The sequence is defined on the set of integers by the recursion formula Xn+1 = Axn+C (mod M), where x, is the nth member of the sequence, and A, C, and M are parameters that can be adjusted for convenience and to ensure the pseudo-random nature of the sequence. For example, M, the modulus, is frequently taken to be the word size on the computer, and A, the multiplier, is chosen to yield both a long period for the sequence and good statistical properties. When M is a power of 2, it has been shown that a suitable sequence can be generated if, among other things, C is odd and A satisfies A = 5 (mod 8) (that is, A-5 is a multiple of 8). A simple example of the generation of a 5-bit number sequence using these conditions would be to set M = 32 (5 bits), A = 21, C = 1, and x0 = 13. This yields the sequence 13, 18, 27, 24, 25, 14, 7, 20, 5, 10..., or, in binary, 01101, 10010, 11011, 11000, 11001, 01110, 00111, 10100, 00101, 01010,.... (1) |