Generating Uniform Random Numbers
This example simulates rolling three dice 10,000 times and plots the distribution of the total:
h = HISTOGRAM(d1 + d2 + d3, LOCATIONS=hlocs)
p = BARPLOT(hlocs, h)
In the above statement, the expression RANDOMU(Seed, 10000) is a 10,000-element, floating-point array of random numbers greater than or equal to 0 and less than 1. Multiplying this array by 6 converts the range to 0 ≤Y < 6. Applying the FIX function yields a 10,000-point integer vector with values from 0 to 5, one less than the numbers on one die. This computation is done three times, once for each die, then the results are added to obtain a vector from 3 to 18.
The HISTOGRAM function makes a vector in which each element contains the number of occurrences of dice rolls whose total is equal to the subscript of the element. Finally, this vector is plotted by the BARPLOT function.
For further information see Additional Examples.
Returns an array of uniformly distributed random numbers of the specified dimensions.
A variable or constant used to initialize the random sequence on input, and in which the state of the random number generator is saved on output.
If Seed is:
- an undefined variable — the generic state is used and the resulting generic state array is returned in Seed. The generic state is initialized on IDL startup using the time-of-day.
- a constant or expression containing a scalar long integer — the value is used to re-initialize the generic state.
- an expression containing an array of long integers of any length — the value is used to re-initialize the generic state.
- a named variable containing a scalar long integer — the scalar value is used to start a new sequence and the resulting state array is returned in Seed.
- a named variable containing an array of long integers of any length — the array is used to start a new sequence and the resulting state array is returned in Seed.
- a named variable that contains the state array from a previous call — it is used to continue the pseudo-random sequence, and the resulting state array is returned in Seed.
For Seed values equal to a non-whole number, such as float, RANDOMU only recognizes the whole number portion of the constant Seed. For example, using 7.0 as the Seed produces the same result as using 7.5 or 7.9.
Note: RANDOMN and RANDOMU use the same sequence, so starting or restarting the sequence for one starts or restarts the sequence for the other. Some IDL routines such as CLUST_WTS use the random number generator, so using them will initialize the seed sequence.
Note: Do not alter the seed array returned by this function. The only valid use for the seed argument is to pass it back to a subsequent call. Changing the seed values will corrupt the random sequence.
Note: Each independent random number sequence should maintain its own state variable. To continue a particular sequence over repeated calls to a procedure, the seed variable should be passed in as an argument to your procedure (or stored in a COMMON block for example).
Either an array or a series of scalar expressions specifying the dimensions of the result. If a single argument is specified, it can be either a scalar expression or an array of up to eight elements. If multiple arguments are specified, they must all be scalar expressions. Up to eight dimensions can be specified. If no dimensions are specified, RANDOMU returns a scalar result.
The formulas for the binomial, gamma, and Poisson distributions are from Section 7.3 of Numerical Recipes in C: The Art of Scientific Computing (Second Edition), Cambridge University Press, 1992.
Set this keyword to a 2-element array, [n, p], to generate random deviates from a binomial distribution. If an event occurs with probability p, with n trials, then the number of times it occurs has a binomial distribution.
Note: For n > 1.0 x 107, you should set the DOUBLE keyword.
Set this keyword to return a double-precision random number.
RANDOMU constructs double-precision uniform random numbers with 53-bits of precision by combining two single-precision numbers:
X = [(A >> 5)*226 + (B >> 6)]*2-53
where A and B are 32-bit integer random numbers in the range 0...4294967295, and the shift-right operator >> is used to remove bits of precision. The resulting values will be 0 ≤ X < 1.
Set this keyword to an integer order i > 0 to generate random deviates from a gamma distribution. The gamma distribution is the waiting time to the ith event in a Poisson random process of unit mean. A gamma distribution of order equal to 1 is the same as the exponential distribution.
Note: For GAMMA > 1.0 x 107, you should set the DOUBLE keyword.
Set this keyword to return integer uniform random deviates in the range [0...231 – 1], using the Mersenne Twister algorithm. If LONG is set, all other keywords (except RAN1) are ignored.
Set this keyword to generate random deviates from a normal distribution.
Set this keyword to the mean number of events occurring during a unit of time. The POISSON keyword returns a random deviate drawn from a Poisson distribution with that mean.
Note: For POISSON > 1.0 x 107, you should set the DOUBLE keyword.
Set this keyword to use the older ran1 algorithm, which was the default in IDL 8.1 and earlier.
Note: The RAN1 keyword is only provided for backwards compatibility, for cases where you need to regenerate the same random sequence from an older version of IDL. Because of the limitations in the algorithm's design, it is strongly recommended that you use the default Mersenne Twister algorithm.
Note: See below for a detailed description of this algorithm.
Set this keyword to return unsigned-integer uniform random deviates in the range [0...232 – 1], using the Mersenne Twister algorithm. If ULONG is set, all other keywords are ignored. The ULONG keyword cannot be used with the RAN1 keyword.
Set this keyword to generate random deviates from a uniform distribution. This is the default behavior.
Using the Seed Argument
If this is the first time that RANDOMU has been called and you start the sequence with an undefined variable, then IDL initializes the sequence with the system time:
; Initialize the sequence with an undefined variable,
; and generate one random number
seed = !NULL
randomValue = RANDOMU(seed)
The new state is saved in seed. To generate repeatable experiments, begin the sequence with a particular seed:
; Initialize the sequence and generate 3 random numbers:
seed = 12345
randomValue = RANDOMU(seed, 3)
; Restart the sequence.
seed = 12345
randomValue = RANDOMU(seed, 3)
0.923121 0.333147 0.197888
0.923121 0.333147 0.197888
To continue generating numbers from the same sequence, pass in the seed variable:
0.949422 0.783800 0.983885
Using an Array of Seed Values
To generate numbers using an array of integers as the seed:
seed = [1798157082, 2109670255, 1881608512, 763029868, 1350847629]
; Pass in the seed array as an expression so our seed is unchanged
x1 = RANDOMU(seed[*], 100000)
; Tweak the last seed value slightly and regenerate
x2 = RANDOMU(seed[*], 100000)
; Should get a completely different random sequence...
To obtain a sequence of 1000 exponential (gamma distribution, order 1) deviates, type:
Result = RANDOMN(seed, 1000, GAMMA=1)
The result contains a random series of waiting times for events occurring an average of one per time period.
To obtain a series of 1000 random elapsed times required for the arrival of two events, type:
;Returns a series of 1000 random elapsed times required for the
;arrival of two events.
Result = RANDOMN(seed, 1000, GAMMA=2)
To obtain a 128 x 128 array filled with Poisson deviates, with a mean of 1.5, type:
Result = RANDOMN(seed, 128, 128, POISSON=1.5)
To simulate the count of “heads” obtained when flipping a coin 10 times, type:
Result = RANDOMN(seed, BINOMIAL=[10,.5])
Notes on the RANDOMU algorithm
- The random number generator uses the Mersenne Twister (MT19937) algorithm.
MT has a very long period of 219937 − 1, and has been optimized for use in Monte Carlo simulations.
- The algorithm is uniformly distributed (to 32-bit accuracy) in any space up to 623 dimensions (the k-distribution test).
- The algorithm is not suitable for cryptography. Since the algorithm is based on linear recursion, observing a sufficiently long sequence could potentially allow you to predict all future iterates.
- IDL uses the mt19937ar.c code, with only one modification: The genrand_real2() function, which generates 32-bit (single-precision) floating point numbers, could occasionally return a value of 1.0 when the return value was cast to a single-precision float. To correct this, if the value 1.0 occurs, it is rejected, and the next random number in the sequence is returned (this is repeated until a value < 1.0 is found).
- The initial seed may be either a scalar long integer or an array of long integers of any length. Passing in an array of integers allows you to select a larger "space" for the seed, and avoids the possibility of accidentally picking the same scalar seed for two different sequences.
- Note that when returning a large number of values, a particular value may occur more than once in the returned array. This is the expected behavior.
- ran1 is a multiplicative congruential algorithm, based off of the algorithm given in Section 7.1 of Numerical Recipes in C: The Art of Scientific Computing (Second Edition), published by Cambridge University Press. To remove low-order serial correlations, a Bays-Durham shuffle is added. During the shuffle, instead of outputting the n-th value in the sequence on the n-th call, it is output on a randomized later call, n+32 on average.
- In IDL's ran1 algorithm, if the value 1.0 occurs, it is rejected, and the next random number in the sequence is returned (this is repeated until a value < 1.0 is found).
- Like all software-based random number generators, the RANDOMU algorithm is a compromise between performance and statistical robustness. Because of the Bays-Durham shuffle, the algorithm should not exhibit any low-order serial correlations. Nevertheless, once the number of generated random numbers approaches the order of the period (approximately 108 or 100,000,000), then the numbers will begin to fail statistical tests.
- Because of the algorithm's design, when returning 32-bit integer random numbers (/LONG is set), a particular value will never be repeated until approximately 231 values have been generated.
|8.2.2||Changed to use the Mersenne Twister algorithm, added the ability to use an array of integers as the seed, added the RAN1 and ULONG keywords.|
The Mersenne Twister algorithm was proposed in the following paper:
M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator", ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, January pp.3-30 (1998). See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for details.
This page has no comments yet. Be the first one!