Bielajew A.F. Fundamentals of the Monte Carlo method for neutral and charged particle transport

Подождите немного. Документ загружается.

2.2. DISCRETE RANDOM VARIABLES 21

• c(x) is a monotonically increasing function of x as a result of p(x)alwaysbeingpositive

and the deﬁnition of c(x) in Equation 2.15.

Cumulative probability distribution functions can be related to uniform random numbers to

provide a way for sampling these distributions. We will complete this discussion in Chapter 4.

Cumulative probability distribution functions for multi-dimensional probability distribution

functions are usually deﬁned in terms of the one-dimensional forms of the marginal and

conditional probability distribution functions.

2.2 Discrete random variables

A more complete discussion of probability theory would include some discussion of discrete

random variables. An example would be the results of ﬂipping a coin or a card game. We

will have some small use for this in Chapter 5 and will introduce what we need at that point.

22 CHAPTER 2. ELEMENTARY PROBABILITY THEORY

Bibliography

Problems

1. Which of the following are candidate probability distributions? For those that are not,

explain. For those that are, determine the normalization constant N. Those that are

proper probability distributions, which contain moments that do not exist?

f(x)=N exp(−µx)0≤ x<∞

f(x)=N exp(−µx)0≤ x<Λ/µ

f(x)=N sin(x)0≤ x<π

f(x)=N sin(x)0≤ x<2π

f(x)=

√

0 ≤ x<1

f(x)=

√

1 ≤ x<∞

f(x)=N

+ a

)

3/2

0 ≤ x<∞

2. Verify that Equations 2.4, 2.6, and 2.7 are true probability distributions.

3. Prove Equation 2.11. Simplify in the case that x and y are independent.

24 BIBLIOGRAPHY

Chapter 3

Random Number Generation

Anyone who considers arithmetical methods of producing random digits is, of

course, in a state of sin.

John von Neumann (1951)

The “pseudo” random number generator (RNG) is the “soul” or “heartbeat” of a Monte

Carlo simulation. It is what generates the pseudo-random nature of Monte Carlo simulations

thereby imitating the true stochastic or random nature of particle interactions. Consequently,

much mathematical study has been devoted to RNG’s [Ehr81, Knu81, Jam88]. These three

references are excellent reviews of RNG theory and methods up to about 1987. The following

references contain more modern material [MZT90, MZ91, L

94, Jam94, Knu97]. It must also

be noted that random number generation is an area of active research. Information that

was given last year may be proven to be misleading this year. The best way to stay in tune

is to track the discussions concerning this topic on the web sites of organizations for which

random number generation is critical. A particularly good one is CERN’s site (www.cern.ch).

CERN is the European Laboratory for Particle Physics. Monte Carlo applications are quite

important in particle physics.

Sometimes one hears the opinion that Monte Carlo codes should be connected somehow to

a source of “true” random numbers. Such true random numbers could be produced by the

noise in an electronic circuit or the time intervals between decays of a radioactive substance.

There are two good reasons NOT to do this. Either a piece of hardware producing the

random numbers would have to be interfaced to a computer somehow, or an array contain-

ing enough random numbers would have to stored. Neither would be practical. However,

the most compelling reason for using mathematically-derived pseudo-random numbers is

repeatability—essential for code debugging. When a Monte Carlo code matures, error dis-

covery becomes less frequent, errors become more subtle and are revealed after the simulation

has run for a long time. Replaying the precise sequence of events that leads to the fault is

essential.

26 CHAPTER 3. RANDOM NUMBER GENERATION

We will not endeavor to explain the theory behind random number generation, merely give

some guidelines for good use. The operative phrase to be used when considering RNG’s

is “use extreme caution”. DO USE an RNG that is known to work well and is widely

tested. DO NOT FIDDLE with RNG’s unless you understand thoroughly the underlying

mathematics and have the ability to test the new RNG thoroughly. DO NOT TRUST

RNG’s that come bundled with standard mathematical packages. For example, DEC’s RAN

RNG (a system utility) and IBM’s RANDU (part of the SSP mathematical package) are known

to give strong triplet correlations. This would aﬀect, for example, the “random” seeding of

an isotropic distribution of point sources in a 3-dimensional object. A picture of an artefact

generated by these RNG’s is given in Figure 3.1. This is known as the “spectral” property

of LCRNG’s.

The gathering of random numbers into planes is a well-known artefact of RNG’s. Marsaglia’s

classic paper [Mar68] entitled “Random numbers fall mainly in the planes”, describes how

random numbers gather into (n − 1)-dimensional hyperplanes in n-space. Good RNG’s ei-

ther maximise the number of planes that are constructed to give the illusion of randomness

or practically eliminate this artefact entirely. One must be aware of this behaviour in case

anomalies do occur. In some cases, despite the shortcoming of RNG’s, no anomalies are de-

tected. An example of this is the same data that produced the obvious artefact in Figure 3.1

but displayed with a 10

◦

rotation about the z-axis does not exhibit the artefact. This is

shown in Figure 3.2.

3.1 Linear congruential random number generators

Most computer architectures support 32-bit 2’s-complement integer arithmetic

. The follow-

ing equation describes a linear congruential random number generator (LCRNG) suitable

for machines that employ 2’s-complement integer arithmetic:

n+1

=mod(aX

+ c, 2

) . (3.1)

This LCRNG generates a 32-bit string of random bits X

n+1

from another representation one

step earlier in the cycle, X

. Upon multiplication or addition, the high-order bits (greater

than position 32) are simply lost leaving the low-order bits scrambled in a pseudo-random

In 32-bit 2’s-complement integer arithmetic

00000000000000000000000000000000 = 0

00000000000000000000000000000001 = 1

00000000000000000000000000000010 = 2 ···

01111111111111111111111111111111 = 2

− 1 = 2147483647

10000000000000000000000000000000 = −2

= −2147483648

10000000000000000000000000000001 = −2

+1=−2147483647

10000000000000000000000000000010 = −2

+2=−2147483646 ···

11111111111111111111111111111111 = −1.

3.1. LINEAR CONGRUENTIAL RANDOM NUMBER GENERATORS 27

Marsaglia planes − View 1

Figure 3.1: The gathering of random numbers into two-dimensional planes when a three-

dimensional cube is seeded.

28 CHAPTER 3. RANDOM NUMBER GENERATION

Marsaglia planes − View 2

Figure 3.2: The identical data in Figure 3.1 but rotated by 10

◦

about the z-axis.

3.1. LINEAR CONGRUENTIAL RANDOM NUMBER GENERATORS 29

fashion. In this equation a is a “magic” multiplier and c is an odd number. The operations

of addition and multiplication take place as if X

n+1

, X

, a,andc are all 32-bit integers.

It is dangerous to consider all the bits as individually random. In fact, the higher-order bits

are more random than the lower-order bits. Random bits are more conveniently produced

by techniques we will discuss later.

The multiplier a is a magic number. Although there are guidelines to follow to determine a

good multiplier, optimum ones are determined experimentally. Particularly good examples

are a = 663608941 and a = 69069. The latter has been suggested by Donald Knuth as the

“best” 32-bit multiplier. It is also easy to remember!

When c is an odd number, the cycle length of the of the LCRNG is 2

(about 4 billion,

in eﬀect, creating every integer in the realm of possibility and an artefactually uniform

random number when converted to ﬂoating point numbers. When c is set to zero, the

LCRNG becomes what is known as a multiplicative congruential random number generator

(MCRNG) with a cycle of 2

(about 1 billion) but with faster execution, saving a fetch and

an integer addition. The seed, X

can be any integer in the case of an LCRNG. In the case

of a MCRNG, it is particularly critical to seed the RNG with an odd number. Conventional

practice is to use either a large odd number (a “random” seed) or a large prime number. If

one seeds the MCRNG with an even integer a reduced cycle length will be obtained. The

severity of the truncation is proportional to the number of times the initial seed can be

divided by two.

The conversion of the random integer to uniform on the range 0 to 1 contains its subtleties as

well and is, to some extent, dependent on architecture. The typical mathematical conversion

is:

=1/2+X

. (3.2)

One subtlety to note is that this should be expected to produce a range 0 ≤ r

< 1 because

the integer representation is asymmetric—it has one more negative integer than positive

integer.

One might code the initialization, generation and conversion in the following fashion:

:Initialisation: i = 987654321

:Iteration: i = i * 663608941

r=0.5+i*0.23283064e-09

Some computers will generate one exact ﬂoating-point zero with this algorithm, others may

generate 2! One way out of this situation is to make the conversion factor slightly smaller, say

0.2328306e-09. This makes the range to be approximately 10

−7

≤ r

< 1 −10

−7

and appears

30 CHAPTER 3. RANDOM NUMBER GENERATION

to work on all computers. Another beneﬁt of this is that it avoids an exact zero which can

cause problems in certain sampling algorithms. We will see more on this later. Although

the endpoints of the distribution are truncated, it usually has no eﬀect. Moreover, if you are

depending upon good uniformity at the endpoints of the random number distribution, then

you probably have to resort to special techniques to obtain them.

3.2 Long sequence random number generators

For many practical applications, cycle lengths of one or four billion are just simply inade-

quate. In fact, a modern workstation just calculating random numbers with these RNG’s

would cycle them in just a few minutes!

One approach is to employ longer integers! The identical algorithm may be employed with

64-bit integers producing a sequence length of 2

or about 1.84 × 10

inthecaseofthe

LCRNG or 2

or about 4.61×10

for the MCRNG. A well-studied multiplier for this purpose

is a = 6364136223846793005. This 64-bit RNG is an excellent choice for the emerging 64-bit

architectures (assuming that they support 64-bit 2’s complement integer arithmetic) or it

can be “faked” using 32-bit architectures [Bie86]. However, this latter approach is no longer

recommended since more powerful long-sequence RNG’s have been developed.

A new approach called the “subtract-with-borrow” algorithm was developed by George

Marsaglia and co-workers [MZT90, MZ91]. This algorithm was very attractive. It was

implemented in ﬂoating-point arithmetic and was portable to all machines. Initialization

and restart capabilities are more involved than for LCRNG’s but the long sequence lengths

and speed of execution make it worthwhile.

Tezuka, l’Ecuyer and Couture [TlC93, Tez94, Cl96] proved that Marsaglia’s algorithm is in

fact equivalent to a linear congruential generator but with a very big integer word length

which greatly reduced the “spectral” anomalies. Since LCRNG’s are so well-studied this is

actually a comfort since for a while the “subtract-with-borrow” algorithm were being em-

ployed with little theoretical understanding, something many researchers felt uncomfortable

with.

The spectral property of this new class of generator has been addressed by L¨uscher [L

94] who

used Kolmogorov’s chaos theory to show how the algorithm could be improved by skipping

some random numbers.

A version of L¨uscher’s algorithm [Jam94] will be distributed with the source library routines

associated with this course. However, keep in mind that the mathematics of random number

generation is still in its infancy. If you take up a Monte Carlo project at some time in the

future, make sure you become well-informed as to the “state-of-the-art”.

An example of Marsaglia and Zaman’s “subtract-with-borrow” algorithm is given as follows: