Page images
[graphic][merged small][merged small][merged small][subsumed][merged small][merged small][merged small][merged small][merged small][subsumed][subsumed][merged small][merged small][merged small][merged small][merged small][merged small][graphic][merged small][merged small][merged small][merged small][merged small][merged small][merged small][merged small][merged small][merged small][ocr errors]
[ocr errors][merged small]


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 are just half-angle (y = x/2) trigonometric identities.)

[merged small][merged small][merged small][merged small][merged small][merged small][ocr errors][merged small][ocr errors][ocr errors][merged small][ocr errors][merged small][ocr errors][merged small][ocr errors][merged small][ocr errors][merged small][merged small]

Now t = tgg,

bution function is
both (independently) equidistributed between 0 and 1, and then

} for

and at


Mr. Stan Ulam

Post Office Box 1663
Santa Fe

New Mexico

Dear Stan:

Thanks for your letter of the 19th. and I are looking forward to the trip and visit at Los Alamos this Summer. I need not tell you that Klari I have already received the necessary papers from Carson Mark. out and returned mine yesterday; Klari's will follow today. I filled

€45°, lies between 0 and 1, and its distri-
Hence one may pick pairs of numbers t, s


With best regards from house to house.


The second method may occasionally be better than the first one.
some cases combinations of both may be hest; e.g., form random pairs
sin x1

cos x

with equidistributed between 0° and 300°. The obvious way consists of
using the sin
- COS tables (with interpolation). This is clearly closely
related to the first method. This is an alternative procedure:
= 2 t
1- t2



2) equidistributed between 0° and 180°.
Then the 7 will have to be replaced randomly by 7,
• .3 7.

with y (which is

This can be done by using random digits 0,

I am very glad that preparations for the random numbers work are to begin soon. In this connection, I would like to mention this: Assume that you have several random number distributions, each equidistributed in Assume that you want one with the


(xi), (y'); (z is,

One way to

distribution function (density) 'f(§) d§ : (§')
form it is to form the cumulative distribution function: g(§) = $$ fl§Id}
to invert it (x1 = § 2
with this hik), or some approximant polynomial. This is, as I see, the
हु x = g 1})
and to form '




method that you have in mind.

An alternative, which works if and all values of §) lie in


0, 1, is this: Scan pairs


y'≤ fix') or not.

= x

to whether
in the second case form no

and use or reject xy
In the first case, put
at that step.



t, s

t но step

(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.

Yours, as ever,


John von Neumann

(1 + t2) SEI



[ocr errors]


[ocr errors][merged small]
[ocr errors]

May 21, 1947

[ocr errors][merged small][merged small][merged small][merged small]
[ocr errors][ocr errors][ocr errors]

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 ≤f(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 x' 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



Fig. 4. If two independent sets of random numbers are used, one of which (x') extends uniformly over the range of the distribution function f and the other (y') extends over the domain of f, then an acceptancerejection technique based on whether or not y'≤ f(x) will generate a distribution for (x') whose density is f(x') dx'.

Reject x1 since y1 > f(x1).




Accept x2 since y2 < f(x2).


X 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

reduction technique. Such techniques 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 x1 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

[ocr errors]

very next free flight are accumulated after each collision.

he Monte Carlo method has proven to be a powerful and useful tool. In fact, "solitaire games" now range from the neutron- and photon-transport codes through the evaluation of multi-dimensional integrals, the exploration of the properties of high-temperature plasmas, and into the quantum mechanics of systems too complex for other methods. Quite a handful. ■


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.

[blocks in formation]
[ocr errors][ocr errors]
« PreviousContinue »