Request for Comments
This is an opensource document; for an updated version, see the source code or its rendering on GitHub. You can send comments on this document either on CodeProject or on the GitHub issues page.
Comments on any aspect of this document are welcome, but especially comments on the following:
 Corrections to any method given on this page.
 Requests to provide an implementation of any method given here in other programming languages, in addition to Python.
 If there is enough interest by readers, I may discuss approaches to generate random mazes, graphs, matrices, or paths.
 Suggestions to trim the size of this document, such as by limiting it to the most common and most useful methods for generating random numbers.
Introduction
This page discusses many ways applications can do random number generation and sampling from an underlying random number generator (RNG), often with pseudocode. Those methods include—
 ways to generate uniform random numbers from an underlying RNG (such as the core method,
RNDINT(N)
),  ways to generate randomized content and conditions, such as Boolean conditions, shuffling, and sampling unique items from a list, and
 generating nonuniform random numbers, including weighted choice, the normal distribution, and other probability distributions.
Sample Python code that implements many of the methods in this document is available.
All the random number methods presented on this page—
 assume the existence of an underlying RNG,
 make no assumptions on the underlying RNG’s implementation (e.g., whether that RNG is a deterministic RNG or some other kind), and
 make no assumptions on the statistical quality or predictability of the underlying RNG.
In general, this document does not cover:
 How to choose an underlying RNG for a particular application, including in terms of security, performance, and quality. I have written more on RNG recommendations in another document.
 Techniques that are specific to an application programming interface.
 Techniques that are specific to certain kinds of RNGs. This includes generating sequences of unique integers using specific kinds of deterministic RNGs.
 Seemingly random numbers that are specifically generated using hash functions, including pseudorandom functions (as opposed to RNGs). But if such a number is used to initialize a deterministic RNG (that is, to serve as its “seed”), then that RNG is generally within the scope of this document.
Contents
 Request for Comments
 Introduction
 Contents
 Notation and Definitions
 Uniform Random Numbers
RNDINT
: Random Integers in [0, N]RNDINTRANGE
: Random Integers in [N, M]RNDU01
: Random Numbers in [0, 1]RNDNUMRANGE
: Random Numbers in [X, Y]RNDINTEXC
: Random Integers in [0, N)RNDINTEXCRANGE
: Random Integers in [N, M)RNDU01OneExc
,RNDU01ZeroExc
, andRNDU01ZeroOneExc
: Random Numbers in [0, 1), (0, 1], or (0, 1)RNDNUMEXCRANGE
: Random Numbers in [X, Y) Uniform Random Bits
 Special Programming Environments
 Randomization Techniques
 General NonUniform Distributions
 Specific NonUniform Distributions
 Dice
 Normal (Gaussian) Distribution
 Binomial Distribution
 Poisson Distribution
 Gamma Distribution
 Negative Binomial Distribution
 von Mises Distribution
 Stable Distribution
 Hypergeometric Distribution
 Multivariate Normal (Multinormal) Distribution
 Random Numbers with a Given Positive Sum
 Multinomial Distribution
 Gaussian and Other Copulas
 Other NonUniform Distributions
 Geometric Sampling
 Conclusion
 Notes
 Appendix
 License
Notation and Definitions
 The pseudocode conventions apply to this document.
 Intervals. The following notation is used for intervals:
 [
a
,b
) means “a
or greater, but less thanb
“.  (
a
,b
) means “greater thana
, but less thanb
“.  (
a
,b
] means “greater thana
and less than or equal tob
“.  [
a
,b
] means “a
or greater andb
or less”.
 [
 Norm. The norm of one or more real numbers is the square root of the sum of squares of those numbers, that is,
sqrt(num1 * num1 + num2 * num2 + ... + numN * numN)
.  Random number generator (RNG). Software and/or hardware that seeks to generate independent numbers that seem to occur by chance and that are approximately uniformly distributed^{(1)}.
 Maximum even parts. The number of maximum even parts is the highest integer
p
such thatp
itself and all factors of1/p
between 0 and 1 are representable in a given floatingpoint number format. For example— the 64bit IEEE 754 binary floatingpoint format (e.g., Java
double
) has 2^{53} (9007199254740992) maximum even parts (see “Generating uniform doubles in the unit interval” in thexoroshiro+
remarks page for further discussion),  the 32bit IEEE 754 binary floatingpoint format (e.g., Java
float
) has 2^{24} (16777216) maximum even parts,  the 64bit IEEE 754 decimal floatingpoint format has 10^{16} maximum even parts, and
 the .NET Framework decimal format (
System.Decimal
) would have 2^{96} maximum even parts, but 2^{96} is not representable, so it actually has 2^{95} maximum even parts instead.
 the 64bit IEEE 754 binary floatingpoint format (e.g., Java
Uniform Random Numbers
This section describes how an underlying RNG can be used to generate uniformlydistributed random numbers. Here is an overview of the methods described in this section.
 Random Integers:
RNDINT
,RNDINTEXC
,RNDINTRANGE
,RNDINTRANGEEXC
.  Random Numbers in 01 Bounded Interval:
RNDU01
,RNDU01ZeroExc
,RNDU01OneExc
,RNDU01ZeroOneExc
.  Other Random Numbers:
RNDNUMRANGE
,RNDNUMRANGEEXC
.
One method, RNDINT
, described next, can serve as the basis for the remaining methods.
RNDINT
: Random Integers in [0, N]
In this document, RNDINT(maxInclusive)
is the core method for generating independent uniform random integers from an underlying RNG, which is called RNG()
in this section. The random integer is in the interval [0, maxInclusive
]. This section explains how RNDINT
can be implemented for two kinds of underlying RNGs; however, the definition of RNDINT
is not limited to those kinds. (These two kinds were chosen because they were the most commonly seen in practice.)
 Method 1: If
RNG()
outputs integers in the interval [0, positiveMODULUS
) (examples ofMODULUS
include 1,000,000 and 6), thenRNDINT(maxInclusive)
can be implemented as in the pseudocode below.^{(2)}  Method 2: If
RNG()
outputs fixedprecision floatingpoint numbers in the interval [0, 1), then finds
, wheres
is the number of maximum even parts for the floatingpoint format, and use Method 1 above, whereMODULUS
iss
andRNG()
isfloor(RNG() * s)
instead.
METHOD RndIntHelperNonPowerOfTwo(maxInclusive) cx = floor(maxInclusive / MODULUS) + 1 while true ret = cx * RNG() // NOTE: If this method is implemented using a fixed // precision type, the addition operation below should // check for overflow and should reject the number // if overflow would result. ret = ret + RNDINT(cx  1) if ret <= maxInclusive: return ret end END METHOD METHOD RndIntHelperPowerOfTwo(maxInclusive) // NOTE: Finds the number of bits minus 1 needed // to represent MODULUS. This will be a constant // here, though. modBits = ln(MODULUS)/ln(2) // Calculate the bit count of maxInclusive bitCount = 0 tempnumber = maxInclusive while tempnumber > 0 // NOTE: If the programming language implements // division with two integers by truncating to an // integer, the division can be used as is without // using a "floor" function. tempnumber = floor(tempnumber / 2) bitCount = bitCount + 1 end while true // Build a number with `bitCount` bits tempnumber = 0 while bitCount > 0 wordBits = modBits rngNumber = RNG() if wordBits > bitCount wordBits = bitCount // Truncate number to 'wordBits' bits // NOTE: If the programming language supports a bitwise // AND operator, the mod operation can be implemented // as "rndNumber AND ((1 << wordBits)  1)" rngNumber = rem(rngNumber, (1 << wordBits)) end tempnumber = tempnumber << wordBits // NOTE: In programming languages that // support the OR operator between two // integers, that operator can replace the // plus operator below. tempnumber = tempnumber + rngNumber bitCount = bitCount  wordBits end // Accept the number if allowed if tempnumber <= maxInclusive: return tempnumber end END METHOD METHOD RNDINT(maxInclusive) // maxInclusive must be 0 or greater if maxInclusive < 0: return error if maxInclusive == 0: return 0 // N equals modulus if maxInclusive == MODULUS  1: return RNG() // NOTE: Finds the number of bits minus 1 needed // to represent MODULUS (if it's a power of 2). // This will be a constant here, though. modBits=ln(MODULUS)/ln(2) // NOTE: The following condition checks if MODULUS // is a power of 2. This will be a constant here, though. isPowerOfTwo=floor(modBits) == modBits // Special cases if MODULUS is a power of 2 if isPowerOfTwo if maxInclusive == 1: return rem(RNG(), 2) if maxInclusive == 3 and modBits >= 2: return rem(RNG(), 4) if maxInclusive == 255 and modBits >= 8: return rem(RNG(), 256) if maxInclusive == 65535 and modBits >=16: return rem(RNG(), 65535) end if maxInclusive > MODULUS  1: if isPowerOfTwo return RndIntHelperPowerOfTwo(maxInclusive) else return RndIntHelperNonPowerOfTwo(maxInclusive) end else // NOTE: If the programming language implements // division with two integers by truncating to an // integer, the division can be used as is without // using a "floor" function. nPlusOne = maxInclusive + 1 maxexc = floor((MODULUS  1) / nPlusOne) * nPlusOne while true ret = RNG() if ret < nPlusOne: return ret if ret < maxexc: return rem(ret, nPlusOne) end end END METHOD
Notes:
 To generate a random number that’s either 1 or 1, the following idiom can be used:
(RNDINT(1) * 2  1)
. To generate a random integer that’s divisible by a positive integer (
DIV
), generate the integer with any method (such asRNDINT
), letX
be that integer, then generateX  rem(X, DIV)
ifX >= 0
, orX  (DIV  rem(abs(X), DIV))
otherwise. (Depending on the method, the resulting integer may be out of range, in which case this procedure is to be repeated.) A random 2dimensional point on an NxM grid can be expressed as a single integer as follows:
 To generate a random NxM point
P
, generateP = RNDINT(N * M  1)
(P
is thus in the interval [0,N * M
)). To convert a point
P
to its 2D coordinates, generate[rem(P, N), floor(P / N)]
. (Each coordinate starts at 0.) To convert 2D coordinates
coord
to an NxM point, generateP = coord[1] * N + coord[0]
.
RNDINTRANGE
: Random Integers in [N, M]
The naïve way of generating a random integer in the interval [minInclusive
, maxInclusive
], shown below, works well for unsigned integers and arbitraryprecision integers.
METHOD RNDINTRANGE(minInclusive, maxInclusive) // minInclusive must not be greater than maxInclusive if minInclusive > maxInclusive: return error return minInclusive + RNDINT(maxInclusive  minInclusive) END METHOD
The naïve approach won’t work as well, though, for signed integer formats if the difference between maxInclusive
and minInclusive
exceeds the highest possible integer for the format. For fixedlength signed integer formats ^{(3)}, such random integers can be generated using the following pseudocode. In the pseudocode below, INT_MAX
is the highest possible integer in the integer format.
METHOD RNDINTRANGE(minInclusive, maxInclusive) // minInclusive must not be greater than maxInclusive if minInclusive > maxInclusive: return error if minInclusive == maxInclusive: return minInclusive if minInclusive==0: return RNDINT(maxInclusive) // Difference does not exceed maxInclusive if minInclusive > 0 or minInclusive + INT_MAX >= maxInclusive return minInclusive + RNDINT(maxInclusive  minInclusive) end while true ret = RNDINT(INT_MAX) // NOTE: If the signed integer format uses two'scomplement // form, use the following line: if RNDINT(1) == 0: ret = 1  ret // NOTE: If the signed integer format has positive and negative // zero, as is the case for Java `float` and // `double` and .NET's implementation of `System.Decimal`, // for example, use the following // three lines instead of the preceding line; // here, zero will be rejected at a 50% chance because zero occurs // twice in both forms. // negative = RNDINT(1) == 0 // if negative: ret = 0  ret // if negative and ret == 0: continue if ret >= minInclusive and ret <= maxInclusive: return ret end END METHOD
Examples:
 To simulate rolling an Nsided die (N greater than 1), generate a random number in the interval [1, N] by
RNDINTRANGE(1, N)
. Generating a random integer with one base10 digit is equivalent to generating
RNDINTRANGE(0, 9)
. Generating a random integer with N base10 digits (where N is 2 or greater) is equivalent to generating
RNDINTRANGE(pow(10, N1), pow(10, N)  1)
.
RNDU01
: Random Numbers in [0, 1]
The idiom RNDINT(X) / X
(called RNDU01()
in this document), generates a random number in the interval [0, 1], where X
is the number of maximum even parts.
For fixedprecision binary floatingpoint numbers with fixed exponent range (such as Java’s double
and float
), the following pseudocode for RNDU01()
can be used instead. It’s based on a technique devised by Allen Downey, who found that dividing a random number by a constant usually does not yield all representable binary floatingpoint numbers in the desired range. In the pseudocode below, SIGBITS
is the binary floatingpoint format’s precision (the number of digits the format can represent without loss; e.g., 53 for Java’s double
).
METHOD RNDU01() e=SIGBITS while true if RNDINT(1)==0: e = e  1 else: break end sig = RNDINT((1 << (SIGBITS  1))  1) if sig==0 and RNDINT(1)==0: e = e + 1 sig = sig + (1 << (SIGBITS  1)) // NOTE: This multiplication should result in // a floatingpoint number; if `e` is sufficiently // small, the number will underflow to 0 return sig * pow(2, e) END METHOD
RNDNUMRANGE
: Random Numbers in [X, Y]
The naïve way of generating a random number in the interval [minInclusive
, maxInclusive
], is shown in the following pseudocode, which generally works well only if the number format can’t be negative or that format uses arbitrary precision.
METHOD RNDNUMRANGE(minInclusive, maxInclusive) if minInclusive > maxInclusive: return error return minInclusive + (maxInclusive  minInclusive) * RNDU01() END
For fixedpoint or floatingpoint number formats with fixed precision (such as Java’s double
and float
), the pseudocode above can overflow if the difference between maxInclusive
and minInclusive
exceeds the maximum possible value for the format. For such formats, the following pseudocode for RNDU01()
can be used instead. In the pseudocode below, NUM_MAX
is the highest possible value for the number format. The pseudocode assumes that the highest possible value is positive and the lowest possible value is negative.
METHOD RNDNUMRANGE(minInclusive, maxInclusive) if minInclusive > maxInclusive: return error if minInclusive == maxInclusive: return minInclusive // Difference does not exceed maxInclusive if minInclusive >= 0 or minInclusive + NUM_MAX >= maxInclusive return minInclusive + (maxInclusive  minInclusive) * RNDU01() end while true ret = RNDU01() * NUM_MAX // NOTE: If the number format has positive and negative // zero, as is the case for Java `float` and // `double` and .NET's implementation of `System.Decimal`, // for example, use the following: negative = RNDINT(1) == 0 if negative: ret = 0  ret if negative and ret == 0: continue // NOTE: For fixedprecision fixedpoint numbers implemented // using two's complement numbers (note 1), use the following line // instead of the preceding three lines, where `QUANTUM` is the // smallest representable positive number in the fixedpoint format: // if RNDINT(1) == 0: ret = (0  QUANTUM)  ret if ret >= minInclusive and ret <= maxInclusive: return ret end END
Note: Monte Carlo integration uses randomization to estimate a multidimensional integral. It involves evaluating a function at N random points in the domain. The estimated integral is the volume of the domain times mean of those N points, and the error in the estimate is that volume times the square root of the (biascorrected sample) variance of the N points (see the appendix). Often quasirandom sequences (also known as lowdiscrepancy sequences, such as Sobol and Halton sequences), often together with an RNG, provide the “random” numbers to sample the function more efficiently. Unfortunately, the methods to produce such sequences are too complicated to show here.
RNDINTEXC
: Random Integers in [0, N)
RNDINTEXC(maxExclusive)
, which generates a random number in the interval [0, maxExclusive
), can be implemented as follows^{(4)}:
METHOD RNDINTEXC(maxExclusive) if maxExclusive <= 0: return error return RNDINT(maxExclusive  1) END METHOD
Note: The following are alternative ways of generating a random integer in the interval [**0,
maxExclusive
):
floor(RNDNUMEXCRANGE(0, maxExclusive))
. Generate
N = floor(RNDU01OneExc()*(maxExclusive))
untilN < maxExclusive
. (The loop is needed because otherwise, rounding error due to the nature of certain floatingpoint formats can result inmaxExclusive
being returned in rare cases.^{(5)})These approaches, though, are recommended only if the programming language—
 supports floatingpoint number types and no other number types (an example is JavaScript),
 is a dialect of SQL, or
 doesn’t support an integer type that is big enough to fit the number
maxExclusive  1
.
RNDINTEXCRANGE
: Random Integers in [N, M)
RNDINTEXCRANGE
returns a random integer in the interval [minInclusive
, maxExclusive
). It can be implemented using RNDINTRANGE
, as the following pseudocode demonstrates.
METHOD RNDINTEXCRANGE(minInclusive, maxExclusive) if minInclusive >= maxExclusive: return error // NOTE: For signed integer formats, replace the following line // with "if minInclusive >=0 or minInclusive + INT_MAX >= // maxExclusive", where `INT_MAX` has the same meaning // as the pseudocode for `RNDINTRANGE`. if minInclusive >=0 return RNDINTRANGE(minInclusive, maxExclusive  1) end while true ret = RNDINTRANGE(minInclusive, maxExclusive) if ret < maxExclusive: return ret end END METHOD
RNDU01OneExc
, RNDU01ZeroExc
, and RNDU01ZeroOneExc
: Random Numbers in [0, 1), (0, 1], or (0, 1)
Three methods related to RNDU01()
can be implemented as follows, where
X
is the number of maximum even parts.

RNDU01OneExc()
, interval [0, 1), can be implemented in one of the following ways: Generate
RNDU01()
in a loop until a number other than 1.0 is generated this way. This is preferred. RNDINT(X  1) / X
.RNDINTEXC(X) / X
.
Note that
RNDU01OneExc()
corresponds toMath.random()
in Java and JavaScript.  Generate
RNDU01ZeroExc()
, interval (0, 1], can be implemented in one of the following ways: Generate
RNDU01()
in a loop until a number other than 0.0 is generated this way. This is preferred. (RNDINT(X  1) + 1) / X
.(RNDINTEXC(X) + 1) / X
.1.0  RNDU01OneExc()
(but this is recommended only if the set of numbersRNDU01OneExc()
could return — as opposed to their probability — is evenly distributed).
 Generate
RNDU01ZeroOneExc()
, interval (0, 1), can be implemented in one of the following ways: Generate
RNDU01()
in a loop until a number other than 0.0 or 1.0 is generated this way. This is preferred. (RNDINT(X  2) + 1) / X
.(RNDINTEXC(X  1) + 1) / X
.
 Generate
RNDNUMEXCRANGE
: Random Numbers in [X, Y)
RNDNUMEXCRANGE
returns a random number in the interval [minInclusive
, maxExclusive
).
It can be implemented using RNDNUMRANGE
, as the following pseudocode demonstrates.
METHOD RNDNUMEXCRANGE(minInclusive, maxExclusive) if minInclusive >= maxExclusive: return error while true ret = RNDNUMRANGE(minInclusive, maxExclusive) if ret < maxExclusive: return ret end END METHOD
Uniform Random Bits
The idiom RNDINT((1 << b)  1)
is a naïve way of generating a uniform random N
bit integer (with maximum 2^{b – 1}).
In practice, memory is usually divided into bytes, or 8bit unsigned integers in the interval [0, 255]. In this case, a byte array or a block of memory can be filled with random bits by setting each byte to RNDINT(255)
.
Special Programming Environments
In certain programming environments it’s often impractical to implement the uniform random number generation methods just described without recurring to other programming languages. These include the following:
 Microsoft Windows batch files (newer versions of which, at least, include a
%RANDOM%
variable which returns a random integer in the interval [0, 65535]). bash
and other shell scripts (some of which include a$RANDOM
variable which returns a random integer in the interval [0, 65535]).
SQL dialects, such as—
 MySQL (which has a
RAND()
akin toRNDU01OneExc()
),  TSQL (which also has a
RAND()
akin toRNDU01OneExc()
),  PL/SQL (which often has
DBMS_RANDOM.VALUE
akin to eitherRNDU01OneExc()
orRNDNUMEXCRANGE
),  PostgreSQL (which has
RANDOM()
), and  SQLite (which sometimes has
RANDOM()
which returns a random integer in the interval [**2^{63}, 2^{63})),
especially within a single query.
 MySQL (which has a
Readers aware of how these environments can support those uniform random number methods should send me a comment.
Moreover, in functional programming languages such as Haskell, it may be important to separate code that directly uses RNGs from other code, usually by rewriting certain functions to take one or more pregenerated random numbers rather than directly using RNDINT
, RNDNUMRANGE
, or another random number generation method presented earlier in this document.
Randomization Techniques
This section describes commonly used randomization techniques, such as shuffling, selection of several unique items, and creating random strings of text.
Boolean Conditions
To generate a condition that is true at the specified probabilities, use
the following idioms in an if
condition:
 True or false with equal probability:
RNDINT(1) == 0
.  True with X percent probability:
RNDINTEXC(100) < X
.  True with probability X/Y:
RNDINTEXC(Y) < X
.  True with odds of X to Y:
RNDINTEXC(X + Y) < X
.  True with probability X, where X is from 0 through 1 (a Bernoulli trial):
RNDU01OneExc() < X
.
Examples:
 True with probability 3/8:
RNDINTEXC(8) < 3
. True with odds of 100 to 1:
RNDINTEXC(101) < 1
. True with 20% probability:
RNDINTEXC(100) < 20
.
Shuffling
The Fisher–Yates shuffle method shuffles a list (puts its items in a random order) such that all permutations (arrangements) of that list are equally likely to occur, assuming the RNG it uses can choose any one of those permutations. However, that method is also easy to write incorrectly (see also Jeff Atwood, “The danger of naïveté“). The following pseudocode is designed to shuffle a list’s contents.
METHOD Shuffle(list) // NOTE: Check size of the list early to prevent // `i` from being less than 0 if the list's size is 0 and // `i` is implemented using an unsigned type available // in certain programming languages. if size(list) >= 2 // Set i to the last item's index i = size(list)  1 while i > 0 // Choose an item ranging from the first item // up to the item given in `i`. Note that the item // at i+1 is excluded. k = RNDINTEXC(i + 1) // The following is wrong since it introduces biases: // "k = RNDINTEXC(size(list))" // The following is wrong since the algorithm won't // choose from among all possible permutations: // "k = RNDINTEXC(i)" // Swap item at index i with item at index k; // in this case, i and k may be the same tmp = list[i] list[i] = list[k] list[k] = tmp // Move i so it points to the previous item i = i  1 end end // NOTE: An implementation can return the // shuffled list, as is done here, but this is not required. return list END METHOD
An important consideration with respect to shuffling is the nature of the underlying RNG, as I discuss in further detail in my RNG recommendation document on shuffling.^{(6)}
Note: In simulation testing, shuffling is used to relabel items from a dataset at random, where each item in the dataset is assigned one of several labels. In such testing—
 one or more statistics that involve the specific labeling of the original dataset’s groups is calculated (such as the difference, maximum, or minimum of means or variances between groups), then
 multiple simulated datasets are generated, where each dataset is generated by—
 merging the groups,
 shuffling the merged dataset, and
 relabeling each item in order such that the number of items in each group for the simulated dataset is the same as for the original dataset, then
 for each simulated dataset, the same statistics are calculated as for the original dataset, then
 the statistics for the simulated datasets are compared with those of the original.
Sampling With Replacement: Choosing a Random Item from a List
To choose a random item from a list—
 whose size is known in advance, use the idiom
list[RNDINTEXC(size(list))]
; or  whose size is not known in advance, generate
RandomKItemsFromFile(file, 1)
, in pseudocode given in a later section (the result will be a 1item list or be an empty list if there are no items).
Choosing an item this way is also known as sampling with replacement.
Notes:
 Generating a random number in the interval [
mn
,mx
) in increments equal tostep
is equivalent to—
 generating a list of all numbers in the interval [
mn
,mx
) of the formmn + step * x
, wherex >= 0
is an integer, then choosing a random item from the list generated this way.
 Bootstrapping is a method of creating a simulated dataset by choosing random items with replacement from an existing dataset until both datasets have the same size. (The simulated dataset can contain duplicates this way.) Usually, multiple simulated datasets are generated this way, one or more statistics, such as the mean, are calculated for each simulated dataset as well as the original dataset, and the statistics for the simulated datasets are compared with those of the original.
Example: Random Character Strings
To generate a random string of characters:
 Generate a list of the letters, digits, and/or other characters the string can have. Examples are given later in this section.

Build a new string whose characters are chosen from that character list. The pseudocode below demonstrates this by creating a list, rather than a string, where the random characters will be held. It also takes the number of characters as a parameter named
size
. (How to convert this list to a text string depends on the programming language and is outside the scope of this page.)METHOD RandomString(characterList, stringSize) i = 0 newString = NewList() while i < stringSize // Choose a character from the list randomChar = characterList[RNDINTEXC(size(characterList))] // Add the character to the string AddItem(newString, randomChar) i = i + 1 end return newString END METHOD
Notes:
 If the list of characters is fixed, the list can be statically created at runtime or compile time, or a string type as provided in the programming language can be used to store the list as a string.
 Instead of individual characters, the list can consist of strings of one or more characters each (e.g., words or syllables), or indeed any other items. (In that case, the individual strings or items should not be stored as a single string).
 Often applications need to generate a string of characters that’s not only random, but also unique. The best way to ensure uniqueness in this case is to store a list (such as a hash table) of strings already generated and to check newly generated strings against that list. Random number generators alone should not be relied on to deliver unique results. If the strings identify database records, file system paths, or other shared resources, special considerations apply, including the need to synchronize access, but are not discussed further in this document.
 Generating “pronounceable” words or words similar to naturallanguage words is generally more sophisticated than shown above. Often, doing so involves Markov chains. A Markov chain models one or more states (for example, individual letters or syllables), and stores the probabilities to transition between these states (e.g., “b” to “e” with a probability of 0.2, or “b” to “b” with a probability of 0.01). A Markov chain modeling random word generation, for example, can include “start” and “stop” states for the start and end of the word, respectively.
Examples of character lists:
 For an alphanumeric string, or string of letters and digits, the characters can be the basic digits “0” to “9” (U+0030U+0039, nos. 4857), the basic upper case letters “A” to “Z” (U+0041U+005A, nos. 6590), and the basic lower case letters “a” to “z” (U+0061U+007A, nos. 96122), as given in the Unicode Standard.
 For a base10 digit string, the characters can be the basic digits only.
 For a base16 digit (hexadecimal) string, the characters can be the basic digits as well as the basic letters “A” to “F” or “a” to “f”.
Sampling Without Replacement: Choosing Several Unique Items
There are several techniques for choosing k
unique items or values uniformly at random from among n
available items or values, depending on such things as whether n
is known and how big n
and k
are.
 If
n
is not known in advance: Use the reservoir sampling method; see theRandomKItemsFromFile
method in the pseudocode below. Although the pseudocode refers to files and lines, the technique applies to any situation when items are retrieved one at a time from a dataset or list whose size is not known in advance.  If items are to be chosen in order:
 If
n
is relatively small, then theRandomKItemsInOrder
method, in the pseudocode below, demonstrates a solution (based on a technique presented in Devroye 1986, p. 620).  If
n
is relatively large, see the item “Ifn
is relatively large”, later.
 If
 If
n
is relatively small (for example, if there are 200 available items, or there is a range of numbers from 0 to 200 to choose from): Do one of the following: Store all the items in a list, shuffle that list, then choose the first
k
items from that list.  If the items are already stored in a list and the list’s order can be changed, then shuffle that list and choose the first
k
items from the shuffled list.  If the items are already stored in a list and the list’s order can’t be changed, then store the indices to those items in another list, shuffle the latter list, then choose the first
k
indices (or the items corresponding to those indices) from the latter list.
 Store all the items in a list, shuffle that list, then choose the first
 If
k
is much smaller thann
and the items are stored in a list whose order can be changed: Do a partial shuffle of that list, then choose the lastk
items from that list. A partial shuffle proceeds as given in the section “Shuffling“, except the partial shuffle stops afterk
swaps have been made (where swapping one item with itself counts as a swap).  If
k
is much smaller thann
andn
is not very large (for example, less than 5000): Do one of the following: Store all the items in a list, do a partial shuffle of that list, then choose the last
k
items from that list.  If the items are already stored in a list and the list’s order can’t be changed, then store the indices to those items in another list, do a partial shuffle of the latter list, then choose the last
k
indices (or the items corresponding to those indices) from the latter list.
 Store all the items in a list, do a partial shuffle of that list, then choose the last
 If
n  k
is much smaller thann
, and the order in which the items are sampled need not be random: If the items are stored in a list whose order can be changed, then proceed as in step 4, except the partial shuffle involves
n  k
swaps and the firstk
items are chosen rather than the lastk
.  Otherwise, if
n
is not very large, then proceed as in step 5, except the partial shuffle involvesn  k
swaps and the firstk
items or indices are chosen rather than the lastk
.
 If the items are stored in a list whose order can be changed, then proceed as in step 4, except the partial shuffle involves
 Otherwise (for example, if 32bit or larger integers will be chosen so that
n
is 2^{32} or isn
is otherwise very large): Create a hash table storing the indices to items already chosen. When a new index to an item is randomly chosen, check the hash table to see if it’s there already. If it’s not there already, add it to the hash table. Otherwise, choose a new random index. Repeat this process untilk
indices were added to the hash table this way. If the items are to be chosen in order, then a red–black tree, rather than a hash table, can be used to store the indices this way; afterk
indices are added to the tree, the indices (and the items corresponding to them) can be retrieved in sorted order. Performance considerations involved in storing data in hash tables or red–black trees, and in retrieving data from them, are outside the scope of this document.
Choosing several unique items as just described is also known as sampling without replacement.
The following pseudocode implements the RandomKItemsFromFile
and RandomKItemsInOrder
methods referred to in this section.
METHOD RandomKItemsFromFile(file, k) list = NewList() j = 0 endOfFile = false while j < k // Get the next line from the file item = GetNextLine(file) // The end of the file was reached, break if item == nothing: endOfFile = true break end AddItem(list, item) j = j + 1 end i = 1 + k while endOfFile == false // Get the next line from the file item = GetNextLine(file) // The end of the file was reached, break if item == nothing: break j = RNDINTEXC(i) if j < k: list[j] = item i = i + 1 end // We shuffle at the end in case k or fewer // lines were in the file, since in that // case the items would appear in the same // order as they appeared in the file // if the list weren't shuffled. This line // can be removed, however, if the items // in the returned list need not appear // in random order. Shuffle(list) return list end METHOD RandomKItemsInOrder(list, k) i = 0 kk = k ret = NewList() n = size(list) while i < n and size(ret) < k u = RNDINTEXC(n  i) if u <= kk AddItem(ret, list[i]) kk = kk  1 end i = i + 1 end return ret END METHOD
Note: Removing
k
random items from a list ofn
items (list
) is equivalent to generating a new
list byRandomKItemsInOrder(list, n  k)
.
Choosing a Random Date/Time
Choosing a random date/time at or between two others is equivalent to—
 converting the two input date/times to an integer or number (here called
date1
anddate2
, wheredate1
represents the earlier date/time anddate2
the other) at the required granularity, for instance, month, day, or hour granularity (the details of such conversion depend on the date/time format and are outside the scope of this document),  generating
newDate = RNDINTRANGE(date1, date2)
ornewDate = RNDNUMRANGE(date1, date2)
, respectively, and  converting
newDate
to a date/time.
If either input date/time was generated as the random date, but that is not desired, the process just given can be repeated until such a date/time is not generated this way.
Generating Random Numbers in Sorted Order
The following pseudocode describes a method that generates random numbers in the interval [0, 1] in sorted order. count
is the number of random numbers to generate this way. The method is based on an algorithm from Bentley and Saxe 1979.
METHOD SortedRandom(count) list = NewList() k = count c = 1.0 while k > 0 c = pow(RNDU01(), 1.0 / k) * c AddItem(list, c) end return list END METHOD
Alternatively, random numbers can be generated (using any method and where the numbers have any distribution and range) and stored in a list, and the list then sorted using a sorting algorithm. Details on sorting algorithms, however, are beyond the scope of this document.
Rejection Sampling
Rejection sampling is a simple and flexible approach for generating random content that
meets certain requirements. To implement rejection sampling:
 Generate the random content (such as a random number) by any method and with any distribution and range.
 If the content doesn’t meet predetermined criteria, go to step 1.
Example criteria include checking any one or more of—
 whether a random number is prime,
 whether a random number is divisible or not by certain numbers,
 whether a random number is not among recently chosen random numbers,
 whether a random number was not already chosen (with the aid of a hash table, redblack tree, or similar structure),
 whether a random point is sufficiently distant from previous random points (with the aid of a KDtree or similar structure),
 whether a random string matches a regular expression, and
 whether a random number is not included in a “blacklist” of numbers.
(KDtrees, hash tables, redblack trees, primenumber testing algorithms, and regular expressions are outside the scope of this document.)
Random Walks
A random walk is a process with random behavior over time. A simple form of random walk involves generating a random number that changes the state of the walk. The pseudocode below generates a random walk of n random numbers, where STATEJUMP()
is the next number to add to the current state (see examples later in this section).
METHOD RandomWalk(n) // Create a new list with an initial state list=[0] // Add 'n' new numbers to the list. for i in 0...n: AddItem(list, list[i] + STATEJUMP()) return list END METHOD
Examples:
 If
STATEJUMP()
isRNDINT(1) * 2  1
, the random walk generates numbers that differ by 1 or 1, chosen at random. If
STATEJUMP()
isRNDNUMRANGE(1, 1)
, the random state is advanced by a random real number in the interval [1, 1]. If
STATEJUMP()
isBinomial(1, p)
, the random walk models a binomial process, where the state is advanced with probabilityp
. If
STATEJUMP()
isBinomial(1, p) * 2  1
, the random walk generates numbers that each differ from the last by 1 or 1 depending on the probabilityp
.Notes:
 A random process can also be simulated by creating a list of random numbers generated the same way. Such a process generally models behavior over time that does not depend on the time or the current state. Examples of this include
Normal(0, 1)
(for modeling Gaussian white noise) andBinomial(1, p)
(for modeling a Bernoulli process, where each number is 0 or 1 depending on the probabilityp
). Some random walks model random behavior at every moment, not just at discrete times. One example is a Wiener process, with random states and jumps that are normally distributed (a process of this kind is also known as Brownian motion). (For a random walk that follows a Wiener process,
STATEJUMP()
isNormal(mu * timediff, sigma * sqrt(timediff))
, wheremu
is the average value per time unit,sigma
is the volatility, andtimediff
is the time difference between samples.) Some random walks model state changes happening at random times. One example is a Poisson process, in which the time between each event is a random exponential variable (the random time is
ln(RNDU01ZeroOneExc()) / rate
, whererate
is the average number of events per time unit; an inhomogeneous Poisson process results ifrate
can vary with the “timestamp” before each event jump.)
General NonUniform Distributions
Some applications need to choose random items or numbers such that some of them are more likely to be chosen than others (a nonuniform distribution). Most of the techniques in this section show how to use the uniform random number methods to generate such random items or numbers.
Weighted Choice
The weighted choice method generates a random item from among a collection of them with separate probabilities of each item being chosen.
The following pseudocode takes a single list weights
, and returns the index of a weight from that list. The greater the weight, the more likely its index will be chosen. (Note that there are two possible ways to generate the random number depending on whether the weights are all integers or can be fractional numbers.) Each weight should be 0 or greater.
METHOD WeightedChoice(weights) if size(weights) == 0: return error sum = 0 // Get the sum of all weights i = 0 while i < size(weights) sum = sum + weights[i] i = i + 1 end // Choose a random integer/number from 0 to less than // the sum of weights. value = RNDINTEXC(sum) // NOTE: If the weights can be fractional numbers, // use this instead: // value = RNDNUMEXCRANGE(0, sum) // Choose the object according to the given value i = 0 lastItem = size(weights)  1 runningValue = 0 while i < size(weights) if weights[i] > 0 newValue = runningValue + weights[i] // NOTE: Includes start, excludes end if value < newValue: return i runningValue = newValue lastItem = i end i = i + 1 end // Last resort (might happen because rounding // error happened somehow) return lastItem END METHOD
Examples:
Assume we have the following list:
["apples", "oranges", "bananas", "grapes"]
, andweights
is the following:[3, 15, 1, 2]
. The weight for “apples” is 3, and the weight for “oranges” is 15. Since “oranges” has a higher weight than “apples”, the index for “oranges” (1) is more likely to be chosen than the index for “apples” (0) with theWeightedChoice
method. The following pseudocode implements how to get a randomly chosen item from the list with that method.index = WeightedChoice(weights) // Get the actual item item = list[index]Assume the weights from example 1 are used and the list contains ranges of numbers instead of strings:
[[0, 5], [5, 10], [10, 11], [11, 13]]
. If a random range is chosen, a random number can be chosen from that range using code like the following:number = RNDNUMEXCRANGE(item[0], item[1])
. (See also “Mixtures of Distributions“.) Assume the weights from example 1 are used and the list contains the following:
[0, 5, 10, 11, 13]
(one more item than the weights). This expresses four ranges, the same as in example 2. After a random index is chosen withindex = WeightedChoice(weights)
, a random number can be chosen from the corresponding range using code like the following:number = RNDNUMEXCRANGE(list[index], list[index + 1])
. (This is how the C++ library expresses a piecewise constant distribution.)
In all the examples above, the weights sum to 21. However, the weights do not mean that, say, when 21 items are selected, the index for “apples” will be chosen exactly 3 times, or the index for “oranges” exactly 15 times, for example. Each number generated by WeightedChoice
is independent from the others, and each weight indicates only a likelihood that the corresponding index will be chosen rather than the other indices. And this likelihood doesn’t change no matter how many times WeightedChoice
is given the same weights. This is called a weighted choice with replacement, which can be thought of as drawing a ball, then putting it back.
Weighted Choice Without Replacement (Multiple Copies)
To implement weighted choice without replacement (which can be thought of as drawing a ball without putting it back), generate an index by WeightedChoice
, and then decrease the weight for the chosen index by 1. In this way, each weight behaves like the number of “copies” of each item. This technique, though, will only work properly if all the weights are integers 0 or greater. The pseudocode below is an example of this.
// Get the sum of weights // (NOTE: This code assumes that `weights` is // a list that can be modified. If the original weights // are needed for something else, a copy of that // list should be made first, but the copying process // is not shown here. This code also assumes that `list`, // a list of items, was already declared earlier and // has at least as many items as `weights`.) totalWeight = 0 i = 0 while i < size(weights) totalWeight = totalWeight + weights[i] i = i + 1 end // Choose as many items as the sum of weights i = 0 items = NewList() while i < totalWeight index = WeightedChoice(weights) // Decrease weight by 1 to implement selection // without replacement. weights[index] = weights[index]  1 AddItem(items, list[index]) i = i + 1 end
Alternatively, if all the weights are integers 0 or greater and their sum is relatively small, create a list with as many copies of each item as its weight, then shuffle that list. The resulting list will be ordered in a way that corresponds to a weighted random choice without replacement.
Note: Weighted choice without replacement can be useful to some applications (particularly some games) that wish to control which random numbers appear, to make the random outcomes appear fairer to users (e.g., to avoid long streaks of good outcomes or of bad outcomes). When used for this purpose, each item represents a different outcome; e.g., “good” or “bad”, and the lists are replenished once no further items can be chosen. However, this kind of sampling should not be used for this purpose whenever information security (ISO/IEC 27000) is involved, including when predicting future random numbers would give a player or user a significant and unfair advantage.
Weighted Choice Without Replacement (Single Copies)
Weighted choice can also choose items from a list, where each item has a separate probability of being chosen and can be chosen no more than once. In this case, after choosing a random index, set the weight for that index to 0 to keep it from being chosen again. The pseudocode below is an example of this.
// (NOTE: This code assumes that `weights` is // a list that can be modified. If the original weights // are needed for something else, a copy of that // list should be made first, but the copying process // is not shown here. This code also assumes that `list`, // a list of items, was already declared earlier and // has at least as many items as `weights`.) chosenItems = NewList() i = 0 // Choose k items from the list while i < k or i < size(weights) index = WeightedChoice(weights) // Set the weight for the chosen index to 0 // so it won't be chosen again weights[index] = 0 // Add the item at the chosen index AddItem(chosenItems, list[index]) end // `chosenItems` now contains the items chosen
The technique presented here can solve the problem of sorting a list of items such that higherweighted items are more likely to appear first.
Continuous Weighted Choice
The continuous weighted choice method generates a random number that follows a continuous probability distribution (here, a piecewise linear distribution).
The pseudocode below takes two lists as follows:
list
is a list of numbers (which need not be integers). If the numbers are arranged in ascending order, which they should, the first number in this list can be returned exactly, but not the last number.weights
is a list of weights for the given numbers (where each number and its weight have the same index in both lists). The greater a number’s weight, the more likely it is that a number close to that number will be chosen. Each weight should be 0 or greater.
METHOD ContinuousWeightedChoice(list, weights) if size(list) <= 0 or size(weights) < size(list): return error if size(list) == 1: return list[0] // Get the sum of all areas between weights sum = 0 areas = NewList() i = 0 while i < size(list)  1 weightArea = abs((weights[i] + weights[i + 1]) * 0.5 * (list[i + 1]  list[i])) AddItem(areas, weightArea) sum = sum + weightArea i = i + 1 end // Choose a random number value = RNDNUMEXCRANGE(0, sum) // Interpolate a number according to the given value i=0 // Get the number corresponding to the random number runningValue = 0 while i < size(list)  1 area = areas[i] if area > 0 newValue = runningValue + area // NOTE: Includes start, excludes end if value < newValue // NOTE: The following line can also read // "interp = RNDU01OneExc()", that is, a new number is generated // within the chosen area rather than using the point // already generated. interp = (value  runningValue) / (newValue  runningValue) retValue = list[i] + (list[i + 1]  list[i]) * interp return retValue end runningValue = newValue end i = i + 1 end // Last resort (might happen because rounding // error happened somehow) return list[size(list)  1] END METHOD
Example: Assume
list
is the following:[0, 1, 2, 2.5, 3]
, andweights
is the following:[0.2, 0.8, 0.5, 0.3, 0.1]
. The weight for 2 is 0.5, and that for 2.5 is 0.3. Since 2 has a higher weight than 2.5, numbers near 2 are more likely to be chosen than numbers near 2.5 with theContinuousWeightedChoice
method.
Mixtures of Distributions
A mixture consists of two or more probability distributions with separate probabilities of being sampled.
To generate random content from a mixture—
 generate
index = WeightedChoice(weights)
, whereweights
is a list of relative probabilities that each distribution in the mixture will be sampled, then  based on the value of
index
, generate the random content from the corresponding distribution.
Examples:
One mixture consists of two normal distributions with two different means: 1 and 1, but the mean 1 normal will be sampled 80% of the time. The following pseudocode shows how this mixture can be sampled:
index = WeightedChoice([80, 20]) number = 0 // If index 0 was chosen, sample from the mean 1 normal if index==0: number = Normal(1, 1) // Else index 1 was chosen, so sample from the mean 1 normal else: number = Normal(1, 1)A hyperexponential distribution is a mixture of exponential distributions, each one with a separate weight and separate rate. An example is below.
index = WeightedChoice([0.6, 0.3, 0.1]) // Rates of the three exponential distributions rates = [0.3, 0.1, 0.05] // Generate an exponential random number with chosen rate number = ln(RNDU01ZeroOneExc()) / rates[index]Choosing a point uniformly at random from a complex shape (in any number of dimensions) is equivalent to sampling uniformly from a mixture of simpler shapes that make up the complex shape (here, the
weights
list holds the content of each simpler shape). (Content is called area in 2D and volume in 3D.) For example, a simple closed 2D polygon can be triangulated, or decomposed into triangles, and a mixture of those triangles can be sampled.^{(7)}For generating a random integer from multiple nonoverlapping ranges of integers—
 each range has a weight of
(mx  mn) + 1
, wheremn
is that range’s minimum andmx
is its maximum, and the chosen range is sampled by generating
RNDINTRANGE(mn, mx)
, wheremn
is the that range’s minimum andmx
is its maximum.For generating random numbers, that may or may not be integers, from nonoverlapping number ranges, each weight is
mx  mn
instead and the number is sampled byRNDNUMEXCRANGE(mn, mx)
instead.
Random Numbers from a Distribution of Data Points
Generating random numbers (or data points) based on how a list of numbers (or data points) is distributed involves a family of techniques called density estimation, which include histograms, kernel density estimation, and Gaussian mixture models. These techniques seek to model the distribution of data points in a given data set, where areas with more points are more likely to be sampled.
 Histograms are sets of one or more bins, which are generally of equal size. Histograms are mixtures, where each bin’s weight is the number of data points in that bin. After a bin is randomly chosen, a random data point that could fit in that bin is generated (that point need not be an existing data point).
 Gaussian mixture models are also mixtures, in this case, mixtures of one or more Gaussian (normal) distributions.
 Kernel distributions are mixtures of sampling distributions, one for each data point. Estimating a kernel distribution is called kernel density estimation. To sample from a kernel distribution:
 Choose one of the numbers or points in the list at random with replacement.
 Add a randomized “jitter” to the chosen number or point; for example, add a separately generated
Normal(0, sigma)
to the chosen number or each component of the chosen point, wheresigma
is the bandwidth^{(8)}.
This document doesn’t detail how to build a density estimation model. Other references on density estimation include a Wikipedia article on multiplevariable kernel density estimation, and a blog post by M. Kay.
Transformations of Random Numbers
Random numbers can be generated by combining and/or transforming one or more random numbers
and/or discarding some of them.
As an example, “Probability and Games: Damage Rolls” by Red Blob Games includes interactive graphics showing score distributions for lowestof, highestof, dropthelowest, and reroll game mechanics.^{(9)} These and similar distributions can be generalized as follows.
Generate two or more random numbers, each with a separate probability distribution, then:
 Highestof: Choose the highest generated number.
 Dropthelowest: Add all generated numbers except the lowest.
 Rerollthelowest: Add all generated numbers except the lowest, then add a number generated randomly by a separate probability distribution.
 Lowestof: Choose the lowest generated number.
 Dropthehighest: Add all generated numbers except the highest.
 Rerollthehighest: Add all generated numbers except the highest, then add a number generated randomly by a separate probability distribution.
 Sum: Add all generated numbers.
 Mean: Find the mean of all generated numbers (see the appendix).
If the probability distributions are the same, then strategies 1 to 3 make higher numbers more likely, and strategies 4 to 6, lower numbers.
Notes: Variants of strategy 4 — e.g., choosing the second, third, or nthlowest number — are formally called second, third, or nthorder statistics distributions, respectively.
Examples:
 The idiom
min(RNDINTRANGE(1, 6), RNDINTRANGE(1, 6))
takes the lowest of two sixsided die results. Due to this approach, 1 is more likely to occur than 6. The idiom
RNDINTRANGE(1, 6) + RNDINTRANGE(1, 6)
takes the result of two sixsided dice (see also “Dice“). Sampling a Bates distribution involves sampling n random numbers by
RNDNUMRANGE(minimum, maximum)
, then finding the mean of those numbers (see the appendix). A compound Poisson distribution models the sum of n random numbers each generated the same way, where n follows a Poisson distribution (e.g.,
n = Poisson(10)
for an average of 10 numbers). A hypoexponential distribution models the sum of n random numbers following an exponential distribution, each with a separate
lamda
parameter (see “Gamma Distribution“).
Random Numbers from an Arbitrary Distribution
Many probability distributions can be defined in terms of any of the following:
 The cumulative distribution function, or CDF, returns, for each number, the probability for a randomly generated variable to be equal to or less than that number; the probability is in the interval [0, 1].
 The probability density function, or PDF, is, roughly and intuitively, a curve of weights 0 or greater, where for each number, the greater its weight, the more likely a number close to that number is randomly chosen.^{(10)}
If a probability distribution’s PDF is known, one of the following techniques, among others, can be used to generate random numbers that follow that distribution.
 Use the PDF to calculate the weights for a number of sample points (usually regularly spaced). Create one list with the sampled points in ascending order (the
list
) and another list of the same size with the PDF’s values at those points (theweights
). Finally generateContinuousWeightedChoice(list, weights)
to generate a random number bounded by the lowest and highest sampled point. This technique can be used even if the area under the PDF isn’t 1. OR  Use inverse transform sampling. Generate
ICDF(RNDU01ZeroOneExc())
, whereICDF(X)
is the distribution’s inverse cumulative distribution function (inverse CDF, or inverse of the CDF) assuming the area under the PDF is 1. OR 
Use rejection sampling. Choose the lowest and highest random number to generate (
minValue
andmaxValue
, respectively) and find the maximum value of the PDF at or between those points (maxDensity
). The rejection sampling approach is then illustrated with the following pseudocode, wherePDF(X)
is the distribution’s PDF (see also Saucier 2000, p. 39). This technique can be used even if the area under the PDF isn’t 1.METHOD ArbitraryDist(minValue, maxValue, maxDensity) if minValue >= maxValue: return error while True x=RNDNUMEXCRANGE(minValue, maxValue) y=RNDNUMEXCRANGE(0, maxDensity) if y < PDF(x): return x end END METHOD
If both a PDF and a uniform random variable in the interval [0, 1) (randomVariable
) are given, then one of the following techniques can be used to generate a random number that follows that distribution:
 Do the same process as method 1, given earlier, except—
 divide the weights in the
weights
list by the sum of all weights, and  use a modified version of
ContinuousWeightedChoice
that usesrandomVariable
rather than generating a new random number. OR
 divide the weights in the
 Generate
ICDF(randomVariable)
, whereICDF(X)
is the distribution’s inverse CDF (see method 2, given earlier).
If the distribution’s CDF is known, generate ICDF(RNDU01ZeroOneExc())
, where ICDF(X)
is the inverse of that CDF.
Note: Further details on inverse transform sampling or on how to find inverses, as well as lists of PDFs and CDFs, are outside the scope of this page.
Censored and Truncated Distributions
To sample from a censored probability distribution, generate a random number from that distribution and—
 if that number is less than a minimum threshold, use the minimum threshold instead, and/or
 if that number is greater than a maximum threshold, use the maximum threshold instead.
To sample from a truncated probability distribution, generate a random number from that distribution and, if that number is less than a minimum threshold and/or higher than a maximum threshold, repeat this process.
Correlated Random Numbers
According to (Saucier 2000), sec. 3.8, to generate two correlated (dependent) random variables—
 generate two independent and identically distributed random variables
x
andy
(for example, twoNormal(0, 1)
variables or twoRNDU01()
variables), and  calculate
[x, y*sqrt(1  rho * rho) + rho * x]
, whererho
is a correlation coefficient in the interval [1, 1] (ifrho
is 0, the variables are uncorrelated).
Other ways to generate correlated random numbers are explained in the section “Gaussian and Other Copulas“.
Specific NonUniform Distributions
This section contains information on some of the most common nonuniform probability distributions.
Dice
The following method generates a random result of rolling virtual dice.^{(11)} It takes three parameters: the number of dice (dice
), the number of sides in each die (sides
), and a number to add to the result (bonus
) (which can be negative, but the result of the subtraction is 0 if that result is greater).
METHOD DiceRoll(dice, sides, bonus) if dice < 0 or sides < 1: return error if dice == 0: return 0 if sides == 1: return dice ret = 0 if dice > 50 // If there are many dice to roll, // use a faster approach, noting that // the diceroll distribution approaches // a "discrete" normal distribution as the // number of dice increases. mean = dice * (sides + 1) * 0.5 sigma = sqrt(dice * (sides * sides  1) / 12) ret = 1 while ret < dice or ret > dice * sides ret = round(Normal(mean, sigma)) end else i = 0 while i < dice ret = ret + RNDINTRANGE(1, sides) i = i + 1 end end ret = ret + bonus if ret < 0: ret = 0 return ret END METHOD
Examples: The result of rolling—
 four sixsided virtual dice (“4d6”) is
DiceRoll(4,6,0)
, three tensided virtual dice, with 4 added (“3d10 + 4”), is
DiceRoll(3,10,4)
, and two sixsided virtual dice, with 2 subtracted (“2d6 – 2”), is
DiceRoll(2,6,2)
.
Normal (Gaussian) Distribution
The normal distribution (also called the Gaussian distribution) can be implemented using the pseudocode below, which uses the polar method ^{(12)} to generate two normallydistributed random numbers:
mu
(μ) is the mean (average), or where the peak of the distribution’s “bell curve” is.sigma
(σ), the standard deviation, affects how wide the “bell curve” appears. The
probability that a normallydistributed random number will be within one standard deviation from the mean is about 68.3%; within two standard deviations (2 timessigma
), about 95.4%; and within three standard deviations, about 99.7%.
METHOD Normal2(mu, sigma) while true a = RNDU01() b = RNDU01() if a != 0 and RNDINT(1) == 0: a = 0  a if b != 0 and RNDINT(1) == 0: b = 0  b c = a * a + b * b if c != 0 and c <= 1 c = sqrt(2 * ln(c) / c) return [a * mu * c + sigma, b * mu * c + sigma] end end END METHOD
Since Normal2
returns two numbers instead of one, but many applications require only one number at a time, a problem arises on how to return one number while storing the other for later retrieval. Ways to solve this problem are outside the scope of this page, however. The name Normal
will be used in this document to represent a method that returns only one normallydistributed random number rather than two.
Alternatively, or in addition, the following method (implementing a ratioofuniforms technique) can be used to generate normally distributed random numbers.
METHOD Normal(mu, sigma) bmp = sqrt(2.0/exp(1.0)) // about 0.8577638849607068 while true a=RNDU01ZeroExc() b=RNDNUMRANGE(bmp,bmp) if b*b <= a * a * 4 * ln(a) return (b * sigma / a) + mu end end END METHOD
Notes:
 In a standard normal distribution,
mu
= 0 andsigma
= 1. If variance is given, rather than standard deviation, the standard deviation (
sigma
) is the variance’s square root.
Binomial Distribution
A random integer that follows a binomial distribution—
 expresses the number of successes that have happened after a given number of independently performed trials
(expressed astrials
below), where the probability of a success in each trial isp
(which ranges from 0, never, to
1, always, and which can be 0.5, meaning an equal chance of success or failure), and  is also known as Hamming distance, if each trial is treated as a “bit” that’s set to 1 for a success and 0 for a failure, and if
p
is 0.5.
METHOD Binomial(trials, p) if trials < 0: return error if trials == 0: return 0 // Always succeeds if p >= 1.0: return trials // Always fails if p <= 0.0: return 0 i = 0 count = 0 // Suggested by Saucier, R. in "Computer // generation of probability distributions", 2000, p. 49 tp = trials * p if tp > 25 or (tp > 5 and p > 0.1 and p < 0.9) countval = 1 while countval < 0 or countval > trials countval = round(Normal(tp, tp)) end return countval end if p == 0.5 while i < trials if RNDINT(1) == 0 // Success count = count + 1 end i = i + 1 end else while i < trials if RNDU01OneExc() < p // Success count = count + 1 end i = i + 1 end end return count END METHOD
Examples:
 If
p
is 0.5, the binomial distribution models the task “Flip N coins, then count the number of heads.” The idiom
Binomial(N, 0.5) >= C
is true if at least C coins, among N coins flipped, show the successful outcome (for example, heads if heads is the successful outcome). The idiom
Binomial(N, 1/S)
models the task “Roll N Ssided dice, then count the number of dice that show the number S.”
Poisson Distribution
The following method generates a random integer that follows a Poisson distribution and is based on Knuth’s method from 1969. In the method—
mean
is the average number of independent events of a certain kind per fixed unit of time or space (for example, per day, hour, or square kilometer), and can be an integer or a noninteger (the method allowsmean
to be 0 mainly for convenience), and the method’s return value gives a random number of such events within one such unit.
METHOD Poisson(mean) if mean < 0: return error if mean == 0: return 0 p = 1.0 // Suggested by Saucier, R. in "Computer // generation of probability distributions", 2000, p. 49 if mean > 9 p = 1.0 while p < 0: p = round(Normal(mean, mean)) return p end pn = exp(mean) count = 0 while true p = p * RNDU01OneExc() if p <= pn: return count count = count + 1 end END METHOD
Gamma Distribution
The following method generates a random number that follows a gamma distribution and is based on Marsaglia and Tsang’s method from 2000. Usually, the number expresses either—
 the lifetime (in days, hours, or other fixed units) of a random component with an average lifetime of
meanLifetime
, or  a random amount of time (in days, hours, or other fixed units) that passes until as many events as
meanLifetime
happen.
Here, meanLifetime
must be an integer or noninteger greater than 0, and scale
is a scaling parameter that is greater than 0, but usually 1.
METHOD GammaDist(meanLifetime, scale) // Needs to be greater than 0 if meanLifetime <= 0 or scale <= 0: return error // Exponential distribution special case if // `meanLifetime` is 1 (see also Devroye 1986, p. 405) if meanLifetime == 1: return ln(RNDU01ZeroOneExc()) * scale d = meanLifetime v = 0 if meanLifetime < 1: d = d + 1 d = d  (1.0 / 3) // NOTE: 1.0 / 3 must be a fractional number c = 1.0 / sqrt(9 * d) while true x = 0 while true x = Normal(0, 1) v = c * x + 1; v = v * v * v if v > 0: break end u = RNDU01ZeroExc() x2 = x * x if u < 1  (0.0331 * x2 * x2): break if ln(u) < (0.5 * x2) + (d * (1  v + ln(v))): break end ret = d * v if meanLifetime < 1 ret = ret * exp(ln(RNDU01ZeroExc()) / meanLifetime) end return ret * scale end
Distributions based on the gamma distribution:
 3parameter gamma distribution:
pow(GammaDist(a, 1), 1.0 / c) * b
, wherec
is another shape parameter.  4parameter gamma distribution:
pow(GammaDist(a, 1), 1.0 / c) * b + d
, whered
is the minimum value.  Exponential distribution:
GammaDist(1, 1.0 / lamda)
orln(RNDU01ZeroOneExc()) / lamda
, wherelamda
is the inverse scale. Usually,lamda
is the probability that an independent event of a given kind will occur in a given span of time (such as in a given day or year), and the random result is the number of spans of time until that event happens. (This distribution is thus useful for modeling a Poisson process.)1.0 / lamda
is the scale (mean), which is usually the average waiting time between two independent events of the same kind.  Erlang distribution:
GammaDist(n, 1.0 / lamda)
. Expresses a sum ofn
exponential random variables with the givenlamda
parameter.
Negative Binomial Distribution
A random integer that follows a negative binomial distribution expresses the number of failures that have happened after seeing a given number of successes (expressed as successes
below), where the probability of a success in each case is p
(where p <= 0
means never, p >= 1
means always, and p = 0.5
means an equal chance of success or failure).
METHOD NegativeBinomial(successes, p) // Needs to be 0 or greater if successes < 0: return error // No failures if no successes or if always succeeds if successes == 0 or p >= 1.0: return 0 // Always fails (NOTE: infinity can be the maximum possible // integer value if NegativeBinomial is implemented to return // an integer) if p <= 0.0: return infinity // NOTE: If 'successes' can be an integer only, // omit the following three lines: if floor(successes) != successes return Poisson(GammaDist(successes, (1  p) / p)) end count = 0 total = 0 if successes == 1 if p == 0.5 while RNDINT(1) == 0: count = count + 1 return count end // Geometric distribution special case (see Saucier 2000) return floor(ln(RNDU01ZeroExc()) / ln(1.0  p)) end while true if RNDU01OneExc() < p // Success total = total + 1 if total >= successes return count end else // Failure count = count + 1 end end END METHOD
Example: If
p
is 0.5 andsuccesses
is 1, the negative binomial distribution models the task “Flip a coin until you get tails, then count the number of heads.”
von Mises Distribution
In the following method, which generates a random number that follows a von Mises distribution, which describes a distribution of circular angles—
mean
is the mean angle,kappa
is a shape parameter, and the method can return a number within π of that mean.
The algorithm below is the Best–Fisher algorithm from 1979 (as described in Devroye 1986 with errata incorporated).
METHOD VonMises(mean, kappa) if kappa < 0: return error if kappa == 0 return RNDNUMEXCRANGE(meanpi, mean+pi) end r = 1.0 + sqrt(4 * kappa * kappa + 1) rho = (r  sqrt(2 * r)) / (kappa * 2) s = (1 + rho * rho) / (2 * rho) while true u = RNDNUMEXCRANGE(1, 1) v = RNDU01ZeroOneExc() z = cos(pi * u) w = (1 + s*z) / (s + z) y = kappa * (s  w) if y*(2  y)  v >=0 or ln(y / v) + 1  y >= 0 if angle<1: angle=1 if angle>1: angle=1 // NOTE: Inverse cosine replaced here // with `atan2` equivalent angle = atan2(sqrt(1w*w),w) if u < 0: angle = angle return mean + angle end end END METHOD
Stable Distribution
As more and more independent and identically distributed random variables are added
together, their distribution tends to a stable distribution,
which resembles a curve with a single peak, but with generally “fatter” tails than the normal distribution. The pseudocode below uses the Chambers–Mallows–Stuck algorithm. The Stable
method, implemented below, takes two parameters:
alpha
is a stability index in the interval (0, 2].beta
is a skewness in the interval [1, 1]; ifbeta
is 0, the curve is symmetric.
METHOD Stable(alpha, beta) if alpha <=0 or alpha > 2: return error if beta < 1 or beta > 1: return error halfpi = pi * 0.5 unif=RNDNUMEXCRANGE(halfpi, halfpi) while unif==halfpi: unif=RNDNUMEXCRANGE(halfpi, halfpi) // Cauchy special case if alpha == 1 and beta == 0: return tan(unif) expo=ln(RNDU01ZeroExc()) c=cos(unif) if alpha == 1 s=sin(unif) return 2.0*((unif*beta+halfpi)*s/c  beta * ln(halfpi*expo*c/(unif*beta+halfpi)))/pi end z=tan(alpha*halfpi)*beta ug=unif+atan2(z, 1)/alpha cpow=pow(c, 1.0 / alpha) return pow(1.0+z*z, 1.0 / (2*alpha))* (sin(alpha*ug)*cpow)* pow(cos(unifalpha*ug)/expo, (1.0  alpha) / alpha) END METHOD
Extended versions of the stable distribution:
 Fourparameter stable distribution:
Stable(alpha, beta) * sigma + mu
, wheremu
is the mean andsigma
is the scale. Ifalpha
andbeta
are 1, the result is a Landau distribution.  “Type 0” stable distribution:
Stable(alpha, beta) * sigma + (mu  sigma * beta * x)
, wherex
isln(sigma)*2.0/pi
ifalpha
is 1, andtan(pi*0.5*alpha)
otherwise.
Hypergeometric Distribution
The following method generates a random integer that follows a hypergeometric distribution.
When a given number of items are drawn at random without replacement from a collection of items
each labeled either 1
or 0
, the random integer expresses the number of items drawn
this way that are labeled 1
. In the method below, trials
is the number of items
drawn at random, ones
is the number of items labeled 1
in the set, and count
is
the number of items labeled 1
or 0
in that set.
METHOD Hypergeometric(trials, ones, count) if ones < 0 or count < 0 or trials < 0 or ones > count or trials > count return error end if ones == 0: return 0 successes = 0 i = 0 currentCount = count currentOnes = ones while i < trials and currentOnes > 0 if RNDINTEXC(currentCount) < currentOnes currentOnes = currentOnes  1 successes = successes + 1 end currentCount = currentCount  1 i = i + 1 end return successes END METHOD
Example: In a 52card deck of AngloAmerican playing cards, 12 of the cards are face
cards (jacks, queens, or kings). After the deck is shuffled and seven cards are drawn, the number
of face cards drawn this way follows a hypergeometric distribution wheretrials
is 7,ones
is
12, andcount
is 52.
Multivariate Normal (Multinormal) Distribution
The following pseudocode calculates a random point in space that follows a multivariate normal (multinormal) distribution. The method MultivariateNormal
takes the following parameters:
 A list,
mu
(μ), which indicates the means to add to each component of the random point.mu
can benothing
, in which case each component will have a mean of zero.  A list of lists
cov
, that specifies a covariance matrix (Σ, a symmetric positive definite NxN matrix, where N is the number of components of the random point).
METHOD Decompose(matrix) numrows = size(matrix) if size(matrix[0])!=numrows: return error // Does a Cholesky decomposition of a matrix // assuming it's positive definite and invertible ret=NewList() for i in 0...numrows submat = NewList() for j in 0...numrows: AddItem(submat, 0) AddItem(ret, submat) end s1 = sqrt(matrix[0][0]) if s1==0: return ret // For robustness for i in 0...numrows ret[0][i]=matrix[0][i]*1.0/s1 end for i in 0...numrows sum=0.0 for j in 0...i: sum = sum + ret[j][i]*ret[j][i] sq=matrix[i][i]sum if sq<0: sq=0 // For robustness ret[i][i]=math.sqrt(sq) end for j in 0...numrows for i in (j + 1)...numrows // For robustness if ret[j][j]==0: ret[j][i]=0 if ret[j][j]!=0 sum=0 for k in 0...j: sum = sum + ret[k][i]*ret[k][j] ret[j][i]=(matrix[j][i]sum)*1.0/ret[j][j] end end end return ret END METHOD METHOD MultivariateNormal(mu, cov) mulen=size(cov) if mu != nothing mulen = size(mu) if mulen!=size(cov): return error if mulen!=size(cov[0]): return error end // NOTE: If multiple random points will // be generated using the same covariance // matrix, an implementation can consider // precalculating the decomposed matrix // in advance rather than calculating it here. cho=Decompose(cov) i=0 ret=NewList() variables=NewList() for j in 0...mulen: AddItem(variables, Normal(0, 1)) while i<mulen nv=Normal(0,1) sum = 0 if mu == nothing: sum=mu[i] for j in 0...mulen: sum=sum+variables[j]*cho[j][i] AddItem(ret, sum) i=i+1 end return ret end
Examples:
 A binormal distribution (twovariable multinormal distribution) can be sampled using the following idiom:
MultivariateNormal([mu1, mu2], [[s1*s1, s1*s2*rho], [rho*s1*s2, s2*s2]])
, wheremu1
andmu2
are the means of the two random variables,s1
ands2
are their standard deviations, andrho
is a correlation coefficient greater than 1 and less than 1 (0 means no correlation). A logmultinormal distribution can be sampled by generating numbers from a multinormal distribution, then applying
exp(n)
to the resulting numbers, wheren
is each number generated this way. A Beckmann distribution can be sampled by calculating the norm of a binormal random pair (see example 1); that is, calculate
sqrt(x*x+y*y)
, wherex
andy
are the two numbers in the binormal pair.
Random Numbers with a Given Positive Sum
Generating N GammaDist(total, 1)
numbers and dividing them by their sum will result in N random numbers that (approximately) sum to total
(see a Wikipedia article). For example, if total
is 1, the numbers will (approximately) sum to 1. Note that in the exceptional case that all numbers are 0, the process should repeat.
The following pseudocode shows how to generate random integers with a given positive sum. (The algorithm for this was presented in Smith and Tromble, “Sampling Uniformly from the Unit Simplex“, 2004.) In the pseudocode below—
 the method
NonzeroIntegersWithSum
returnsn
positive integers that sum tototal
,  the method
IntegersWithSum
returnsn
nonnegative integers that sum tototal
, and Sort(list)
sorts the items inlist
in ascending order (note that details on sort algorithms are outside the scope of this document).
METHOD NonzeroIntegersWithSum(n, total) if n <= 0 or total <=0: return error ls = NewList() ret = NewList() AddItem(ls, 0) while size(ls) < n c = RNDINTEXCRANGE(1, total) found = false j = 1 while j < size(ls) if ls[j] == c found = true break end j = j + 1 end if found == false: AddItem(ls, c) end Sort(ls) AddItem(ls, total) for i in 1...size(ls): AddItem(ret, list[i]  list[i  1]) return ret END METHOD METHOD IntegersWithSum(n, total) if n <= 0 or total <=0: return error ret = NonzeroIntegersWithSum(n, total + n) for i in 0...size(ret): ret[i] = ret[i]  1 return ret END METHOD
Notes:
 The problem of generating N random numbers with a given positive sum
sum
is equivalent to the problem of generating a uniformly distributed point inside an Ndimensional simplex (simplest convex figure) whose edges all have a length ofsum
units. Generating
N
random numbers with a given positive averageavg
is equivalent to generatingN
random numbers with the sumN * avg
. Generating
N
random numbersmin
or greater and with a given positive sumsum
is equivalent to generatingN
random numbers with the sumsum  n * min
, then addingmin
to each number generated this way. The similar Dirichlet distribution models n random numbers, and can be sampled by generating n+1 random gammadistributed numbers, each with separate parameters, and dividing all those numbers except the last by the sum of all n+1 numbers (see Devroye 1986, p. 594).
Multinomial Distribution
The multinomial distribution models the number of times each of several mutually exclusive events happens among a given number of trials, where each event can have a separate probability of happening. The pseudocode below is of a method that takes two parameters: trials
, which is the number of trials, and weights
, which are the relative probabilities of each event. The method tallies the events as they happen and returns a list (with the same size as weights
) containing the number of successes for each event.
METHOD Multinomial(trials, weights) if trials < 0: return error // create a list of successes list = NewList() for i in 0...size(weights): AddItem(list, 0) for i in 0...trials // Choose an index index = WeightedChoice(weights) // Tally the event at the chosen index list[index] = list[index] + 1 end return list END METHOD
Gaussian and Other Copulas
A copula is a distribution describing the correlation (dependence) between random numbers.
One example is a Gaussian copula; this copula is sampled by sampling from a multinormal distribution, then converting the resulting numbers to uniformlydistributed, but correlated, numbers. In the following pseudocode, which implements a Gaussian copula:
 The parameter
covar
is the covariance matrix for the multinormal distribution. erf(v)
is the error function of the variablev
(see the appendix).
METHOD GaussianCopula(covar) mvn=MultivariateNormal(nothing, covar) for i in 0...size(covar) // Apply the normal distribution's CDF // to get uniform variables mvn[i] = (erf(mvn[i]/(sqrt(2)*sqrt(covar[i][i])))+1)*0.5 end return mvn END METHOD
Each of the resulting uniform numbers will be in the interval [0, 1], and each one can be further transformed to any other probability distribution (which is called a marginal distribution here) by one of the methods given in “Random Numbers from an Arbitrary Distribution“. (See also Cario and Nelson 1997.)
Example: To generate two correlated uniform variables with a Gaussian copula, generate
GaussianCopula([[1, rho], [rho, 1]])
, whererho
is the Pearson correlation coefficient, in the interval [1, 1]. (Note that rank correlation parameters, which can be converted torho
, can better describe the correlation thanrho
itself. For example, for a twovariable Gaussian copula, the Spearman coefficientsrho
can be converted torho
byrho = sin(srho * pi / 6) * 2
. Rank correlation parameters are not further discussed in this document.)
Other kinds of copulas describe different kinds of correlation between random numbers. Examples of other copulas are—
 the Fréchet–Hoeffding upper bound copula [x, x, …, x] (e.g.,
[x, x]
), wherex = RNDU01()
,  the Fréchet–Hoeffding lower bound copula
[x, 1.0  x]
wherex = RNDU01()
,  the product copula, where each number is a separately generated
RNDU01()
(indicating no correlation between the numbers), and  the Archimedean copulas, described by M. Hofert and M. Mächler (2011)^{(13)}.
Other NonUniform Distributions
Most commonly used:
 Beta distribution (
BetaDist(a, b)
):x / (x + GammaDist(b, 1))
, wherex
isGammaDist(a + Poisson(nc * 0.5), 1)
,a
andb
are
the two parameters of the beta distribution, andnc
is a parameter such thatnc
other than 0 indicates a noncentral distribution. The range of the beta distribution is [0, 1).  Cauchy (Lorentz) distribution:
scale * tan(pi * (RNDU01OneExc()0.5)) + mu
, wherescale
is the scale andmu
is the location of the distribution’s curve peak (mode). This distribution is similar to the normal distribution, but with “fatter” tails.  Chisquared distribution:
GammaDist(df * 0.5 + Poisson(sms * 0.5), 2)
, wheredf
is the number of degrees of freedom andsms
is the sum of mean squares (wheresms
other than 0 indicates a noncentral distribution).  Extreme value distribution:
a  ln(ln(RNDU01ZeroOneExc())) * b
, whereb
is the scale anda
is the location of the distribution’s curve peak (mode). This expresses a distribution of maximum values.  Geometric distribution:
NegativeBinomial(1, p)
, wherep
has the same meaning
as in the negative binomial distribution. As used here, this is the number of failures that have happened before a success happens. (Saucier 2000, p. 44, also mentions an alternative definition that includes the success.)  Gumbel distribution:
a + ln(ln(RNDU01ZeroOneExc())) * b
, whereb
is the scale anda
is the location of the distribution’s curve peak (mode). This expresses a distribution of minimum values.  Inverse gamma distribution:
b / GammaDist(a, 1)
, wherea
andb
have the
same meaning as in the gamma distribution. Alternatively,1.0 / (pow(GammaDist(a, 1), 1.0 / c) / b + d)
, wherec
andd
are shape and location parameters, respectively.  Laplace (double exponential) distribution:
(ln(RNDU01ZeroExc())  ln(RNDU01ZeroExc())) * beta + mu
, wherebeta
is the scale andmu
is the mean.  Logarithmic distribution:
min + (max  min) * RNDU01OneExc() * RNDU01OneExc()
, wheremin
is the minimum value andmax
is the maximum value (Saucier 2000, p. 26). In this distribution, numbers closer tomin
are exponentially more likely than numbers closer tomax
.  Logarithmic normal distribution:
exp(Normal(mu, sigma))
, wheremu
andsigma
have the same meaning as in the normal distribution.  Pareto distribution:
pow(RNDU01ZeroOneExc(), 1.0 / alpha) * minimum
, wherealpha
is the shape andminimum
is the minimum.  Rayleigh distribution:
a * sqrt(2 * ln(RNDU01ZeroExc()))
, wherea
is the scale and is greater than 0. Ifa
follows a logarithmic normal distribution, the result is a Suzuki distribution.  Student’s tdistribution:
Normal(cent, 1) / sqrt(GammaDist(df * 0.5, 2 / df))
, wheredf
is the number of degrees of freedom, and cent is the mean of the normallydistributed random number. Acent
other than 0 indicates a noncentral distribution.  Triangular distribution:
ContinuousWeightedChoice([startpt, midpt, endpt], [0, 1, 0])
. The distribution starts atstartpt
, peaks atmidpt
, and ends atendpt
.  Weibull distribution:
b * pow(ln(RNDU01ZeroExc()),1.0 / a) + loc
, wherea
is the shape,b
is the scaleloc
is the location, anda
andb
are greater than 0.
Miscellaneous:
 Arcsine distribution:
min + (max  min) * BetaDist(0.5, 0.5)
, wheremin
is the minimum value andmax
is the maximum value (Saucier 2000, p. 14).  Beta binomial distribution:
Binomial(trials, BetaDist(a, b))
, wherea
andb
are
the two parameters of the beta distribution, andtrials
is a parameter of the binomial distribution.  BetaPERT distribution:
startpt + size * BetaDist(1.0 + (midpt  startpt) * shape / size, 1.0 + (endpt  midpt) * shape / size)
. The distribution starts atstartpt
, peaks atmidpt
, and ends atendpt
,size
isendpt  startpt
, andshape
is a shape parameter that’s 0 or greater, but usually 4. If the mean (mean
) is known rather than the peak,midpt = 3 * mean / 2  (startpt + endpt) / 4
.  Beta prime distribution:
pow(GammaDist(a, 1), 1.0 / alpha) * scale / pow(GammaDist(b, 1), 1.0 / alpha)
, wherea
,b
, andalpha
are shape parameters andscale
is the scale. If a is 1, the result is a Singh–Maddala distribution; if b is 1, a Dagum distribution; if a and b are both 1, a logarithmic logistic distribution.  Beta negative binomial distribution:
NegativeBinomial(successes, BetaDist(a, b))
, wherea
andb
are
the two parameters of the beta distribution, andsuccesses
is a parameter of the negative binomial distribution. If successes is 1, the result is a Waring–Yule distribution.  Birnbaum–Saunders distribution:
pow(sqrt(4+x*x)+x,2)/(4.0*lamda)
, wherex = Normal(0,gamma)
,gamma
is a shape parameter, andlamda
is a scale parameter.  Chi distribution:
sqrt(GammaDist(df * 0.5, 2))
, wheredf
is the number of degrees of freedom.  Cosine distribution:
min + (max  min) * atan2(x, sqrt(1  x * x)) / pi
, wherex = RNDNUMRANGE(1, 1)
andmin
is the minimum value andmax
is the maximum value (Saucier 2000, p. 17; inverse sine replaced withatan2
equivalent).  Double logarithmic distribution:
min + (max  min) * (0.5 + (RNDINT(1) * 2  1) * 0.5 * RNDU01OneExc() * RNDU01OneExc())
, wheremin
is the minimum value andmax
is the maximum value (see also Saucier 2000, p. 15, which shows the wrong X axes).  Fréchet distribution:
b*pow(ln(RNDU01ZeroExc()),1.0/a) + loc
, wherea
is the shape,b
is the scale, andloc
is the location of the distribution’s curve peak (mode). This expresses a distribution of maximum values.  Generalized extreme value (Fisher–Tippett) distribution:
a  (pow(ln(RNDU01ZeroOneExc()), c)  1) * b / c
ifc != 0
, ora  ln(ln(RNDU01ZeroOneExc())) * b
otherwise, whereb
is the scale,a
is the location of the distribution’s curve peak (mode), andc
is a shape parameter. This expresses a distribution of maximum values.  Generalized Tukey lambda distribution:
(s1 * (pow(x, lamda1)1.0)/lamda1  s2 * (pow(1.0x, lamda2)1.0)/lamda2) + loc
, wherex
isRNDU01()
,lamda1
andlamda2
are shape parameters,s1
ands2
are scale parameters, andloc
is a location parameter.  Halfnormal distribution. Parameterizations include:
 Mathematica:
abs(Normal(0, sqrt(pi * 0.5) / invscale)))
, whereinvscale
is a parameter of the halfnormal distribution.  MATLAB:
abs(Normal(mu, sigma)))
, wheremu
andsigma
are the same as in the normal distribution.
 Mathematica:
 Inverse chisquared distribution:
df * scale / (GammaDist(df * 0.5, 2))
, wheredf
is the number of degrees of freedom andscale
is the scale, usually1.0 / df
.  Inverse Gaussian distribution (Wald distribution): Generate
n = mu + (mu*mu*y/(2*lamda))  mu * sqrt(4 * mu * lamda * y + mu * mu * y * y) / (2 * lamda)
, wherey = pow(Normal(0, 1), 2)
, then returnn
ifRNDU01OneExc() <= mu / (mu + n)
, ormu * mu / n
otherwise.mu
is the mean andlamda
is the scale; both parameters are greater than 0. Based on method published in Devroye 1986.  Kumaraswamy distribution:
min + (max  min) * pow(1pow(RNDU01ZeroExc(),1.0/b),1.0/a)
, wherea
andb
are shape parameters,min
is the minimum value, andmax
is the maximum value.  Lévy distribution:
sigma * 0.5 / GammaDist(0.5, 1) + mu
, wheremu
is the location andsigma
is the dispersion.  Logarithmic series distribution:
floor(1.0 + ln(RNDU01ZeroExc()) / ln(1.0  pow(1.0  param, RNDU01ZeroOneExc())))
, whereparam
is a number greater than 0 and less than 1. Based on method described in Devroye 1986.  Logistic distribution:
(ln(x)ln(1.0  x)) * scale + mean
, wherex
isRNDU01ZeroOneExc()
andmean
andscale
are the mean and the scale, respectively.  Maxwell distribution:
scale * sqrt(GammaDist(1.5, 2))
, wherescale
is the scale.  Parabolic distribution:
min + (max  min) * BetaDist(2, 2)
, wheremin
is the minimum value andmax
is the maximum value (Saucier 2000, p. 30).  Pascal distribution:
NegativeBinomial(successes, p) + successes
, wheresuccesses
andp
have the same meaning as in the negative binomial distribution, exceptsuccesses
is always an integer.  Pearson VI distribution:
GammaDist(v, 1) / (GammaDist(w, 1))
, wherev
andw
are shape parameters greater than 0 (Saucier 2000, p. 33; there, an additionalb
parameter is defined, but that parameter is canceled out in the source code).  Power distribution:
pow(RNDU01ZeroOneExc(), 1.0 / alpha) / b
, wherealpha
is the shape andb
is the domain. Nominally in the interval (0, 1).  Power law distribution:
pow(pow(mn,n+1) + (pow(mx,n+1)  pow(mn,n+1)) * RNDU01(), 1.0 / (n+1))
, wheren
is the exponent,mn
is the minimum, andmx
is the maximum. Reference.  Skellam distribution:
Poisson(mean1)  Poisson(mean2)
, wheremean1
andmean2
are the means of the two Poisson variables.  Skewed normal distribution:
Normal(0, x) + mu + alpha * abs(Normal(0, x))
, wherex
issigma / sqrt(alpha * alpha + 1.0)
,mu
andsigma
have the same meaning as in the normal distribution, andalpha
is a shape parameter.  Snedecor’s (Fisher’s) Fdistribution:
GammaDist(m * 0.5, n) / (GammaDist(n * 0.5 + Poisson(sms * 0.5)) * m, 1)
, wherem
andn
are the numbers of degrees of freedom of two random numbers with a chisquared distribution, and ifsms
is other than 0, one of those distributions is noncentral with sum of mean squares equal tosms
.  Tukey lambda distribution:
(pow(x, lamda)pow(1.0x,lamda))/lamda
, wherex
isRNDU01()
andlamda
is a shape parameter (if 0, the result is a logistic distribution).  Zeta distribution: Generate
n = floor(pow(RNDU01ZeroOneExc(), 1.0 / r))
, and ifd / pow(2, r) < (d  1) * RNDU01OneExc() * n / (pow(2, r)  1.0)
, whered = pow((1.0 / n) + 1, r)
, repeat this process. The parameterr
is greater than 0. Based on method described in Devroye 1986. A zeta distribution truncated by rejecting random values greater than some positive integer is called a Zipf distribution or Estoup distribution. (Note that Devroye uses “Zipf distribution” to refer to the untruncated zeta distribution.)
The Python sample code also contains implementations of the power normal distribution, the power lognormal distribution, the negative multinomial distribution, the multivariate tdistribution, the multivariate tcopula, and the multivariate Poisson distribution.
Geometric Sampling
This section contains various geometric sampling techniques.
Random Point Inside a Triangle
The following pseudocode, which
generates, uniformly at random, a point inside a 2dimensional triangle,
takes three parameters, p0
, p1
, and p2
, each of which is a 2item list containing the X and Y
coordinates, respectively, of one vertex of the triangle.
METHOD RandomPointInTriangle( x1=p1[0]p0[0] y1=p1[1]p0[1] x2=p2[0]p0[0] y2=p2[1]p0[1] den=(x1*y2x2*y1) // Avoid division by zero if den==0: den=0.0000001 r=RNDU01() s=RNDU01() xv=r*x1 + s*x2 yv=r*y1 + s*y2 a=(xv*y2  x2*yv)/den b=(x1*yv  xv*y1)/den if a<=0 or b<=0 or a+b>=1 return [x1+x2+p0[0]xv,y1+y2+p0[1]yv] else return [p0[0]+xv, p0[1]+yv] end end
Random Points on the Surface of a Hypersphere
To generate, uniformly at random, an Ndimensional point on the surface of an Ndimensional hypersphere of radius R, generate N Normal(0, 1)
random numbers, then multiply them by R / X
, where X
is those numbers’ norm (if X
is 0, the process should repeat). A hypersphere’s surface is formed by all points lying R units away from a common point in Ndimensional space. Based on a technique described in MathWorld.
This problem is equivalent to generating, uniformly at random, a unit vector (vector with length 1) in Ndimensional space.
Random Points Inside a Ball or Shell
To generate, uniformly at random, an Ndimensional point inside an Ndimensional ball of radius R, either—
 generate N
Normal(0, 1)
random numbers, generateX = sqrt( S  ln(RNDU01ZeroExc()))
, whereS
is the sum of squares of the random numbers, and multiply each random number byR / X
(ifX
is 0, the process should repeat), or  generate N
RNDNUMRANGE(R, R)
random numbers^{(14)} until their norm is R or less,
although the former method “may … be slower” “in practice”, according to a MathWorld article, which was the inspiration for the two methods given here.
To generate, uniformly at random, a point inside an Ndimensional spherical shell (a hollow ball) with inner radius A and outer radius B (where A is less than B), either—
 generate, uniformly at random, a point for a ball of radius B until the norm of the numbers making up that point is A or greater;
 for 2 dimensions, generate, uniformly at random, a point on the surface of a circle with radius equal to
sqrt(RNDNUMRANGE(0, B * B  A * A) + A * A)
(Dupree and Fraley 2004); or  for 3 dimensions, generate, uniformly at random, a point on the surface of a sphere with radius equal to
pow(RNDNUMRANGE(0, pow(B, 3)  pow(A, 3)) + pow(A, 3), 1.0 / 3.0)
(Dupree and Fraley 2004).
Random Latitude and Longitude
To generate, uniformly at random, a point on the surface of a sphere in the form of a latitude and longitude (in radians with west and south coordinates negative)—
 generate the longitude
RNDNUMEXCRANGE(pi, pi)
, where the longitude ranges from π to π, and  generate the latitude
atan2(sqrt(1  x * x), x)  pi / 2
, where—x = RNDNUMRANGE(1, 1)
and the latitude ranges from π/2 to π/2 (the range includes the poles, which have many equivalent forms), orx = 2 * RNDU01ZeroOneExc()  1
and the latitude ranges from π/2 to π/2 (the range excludes the poles).
Reference: “Sphere Point Picking” in MathWorld (replacing inverse cosine with atan2
equivalent).
Conclusion
This page discussed many ways applications can extract random numbers
from random number generators.
I acknowledge the commenters to the CodeProject version of this page, including George Swan, who referred me to the reservoir sampling method.
Notes
^{(1)} An RNG meeting this definition can—
 seek to generate random numbers that are at least costprohibitive (but not necessarily impossible) to predict,
 merely seek to generate number sequences likely to pass statistical tests of randomness,
 be initialized automatically before use,
 be initialized with an applicationspecified “seed”,
 use a deterministic algorithm for random number generation,
 rely, at least primarily, on one or more nondeterministic sources for random number
generation (including by extracting uniformly distributed bits from two or more such sources), or  have two or more of the foregoing properties.
If the software and/or hardware uses a nonuniform distribution, but otherwise meets this definition, it can be converted to use a uniform distribution, at least in theory, using unbiasing or randomness extraction methods that it is outside the scope of this document to describe.
^{(2)} For an exercise solved by this method, see A. Koenig and B. E. Moo, Accelerated C++, 2000; see also a blog post by Johnny Chan.
Note that if MODULUS
is a power of 2 (for example, 256 or 2^{32}), the RNDINT
implementation given may leave unused bits (for example, when truncating a random number to wordBits
bits or in the special cases at the start of the method). How a more sophisticated implementation may save those bits for later reuse is beyond this page’s scope.
^{(3)} This number format describes Bbit signed integers with minimum value 2^{B1} and maximum value 2^{B1} – 1, where B is a positive even number of bits; examples include Java’s short
, int
, and long
, with 16, 32, and 64 bits, respectively. A signed integer is an integer that can be positive, zero, or negative. In two’scomplement form, nonnegative numbers have the highest (most significant) bit set to zero, and negative numbers have that bit (and all bits beyond) set to one, and a negative number is stored in such form by swapping the bits of a number equal to that number’s absolute value minus 1.
^{(4)} RNDINTEXC
is not given as the core random generation method because it’s harder to fill integers in popular integer formats with random bits with this method.
^{(5)} In situations where loops are not possible, such as within an SQL query, the idiom min(floor(RNDU01OneExc() * maxExclusive, maxExclusive  1))
returns an integer in the interval [0, maxExclusive
); however, such an idiom can have a slight, but for most purposes negligible, bias toward maxExclusive  1
.
^{(6)} It suffices to say here that in general, whenever a deterministic RNG is otherwise called for, such an RNG is good enough for shuffling a 52item list if its period is 2^{226} or greater. (The period is the maximum number of values in a generated sequence for a deterministic RNG before that sequence repeats.)
^{(7)} A convex polygon can be trivially decomposed into triangles that have one vertex in common and each have two other adjacent vertices of the original polygon. Triangulation of other polygons is nontrivial and outside the scope of this document.
^{(8)} “Jitter”, as used in this step, follows a distribution formally called a kernel, of which the normal distribution is one example. Bandwidth should be as low or as high as allows the estimated distribution to fit the data and remain smooth. A more complex kind of “jitter” (for multicomponent data points) consists of a point generated from a multinormal distribution with all the means equal to 0 and a covariance matrix that, in this context, serves as a bandwidth matrix. “Jitter” and bandwidth are not further discussed in this document.
^{(9)} That article also mentions a criticalhit distribution, which is actually a mixture of two distributions: one roll of dice and the sum of two rolls of dice.
^{(10)} More formally—
 the PDF is the derivative (instantaneous rate of change) of the distribution’s CDF (that is, PDF(x) = CDF′(x)), and
 the CDF is also defined as the integral of the PDF,
provided the PDF’s values are all 0 or greater and the area under the PDF’s curve is 1.
^{(11)} The “Dice” section used the following sources:
 Red Blob Games, “Probability and Games: Damage Rolls” was the main source for the diceroll distribution. The method
random(N)
in that document corresponds toRNDINTEXC(N)
in this document.  The MathWorld article “Dice” provided the mean of the dice roll distribution.
 S. Eger, “Stirling’s approximation for central extended binomial coefficients”, 2014, helped suggest the variance of the dice roll distribution.
^{(12)} The method that formerly appeared here is the BoxMuller transformation: mu + radius * cos(angle)
and mu + radius * sin(angle)
, where angle = 2 * pi * RNDU01OneExc()
and radius = sqrt(2 * ln(RNDU01ZeroExc())) * sigma
, are two independent normallydistributed random numbers. A method of generating approximate standard normal random numbers, which consists of summing twelve RNDU01OneExc()
numbers and subtracting by 6 (see also “Irwin–Hall distribution” on Wikipedia), results in values not less than 6 or greater than 6; on the other hand, in a standard normal distribution, results less than 6 or greater than 6 will occur only with a generally negligible probability.
^{(13)} Hofert, M., and Maechler, M. “Nested Archimedean Copulas Meet R: The nacopula Package”. Journal of Statistical Software 39(9), 2011, pp. 120.
^{(14)} The N numbers generated this way will form a point inside an Ndimensional hypercube with length 2 * R
in each dimension and centered at the origin of space.
Appendix
Implementation of erf
The pseudocode below shows how the error function erf
can be implemented, in case the programming language used doesn’t include a builtin version of erf
(such as JavaScript at the time of this writing). In the pseudocode, EPSILON
is a very small number to end the iterative calculation.
METHOD erf(v) if v==0: return 0 if v<0: return erf(v) if v==infinity: return 1 // NOTE: For Java `double`, the following // line can be added: // if v>=6: return 1 i=1 ret=0 zp=(v*v) zval=1.0 den=1.0 while i < 100 r=v*zval/den den=den+2 ret=ret+r // NOTE: EPSILON can be pow(10,14), // for example. if abs(r)<EPSILON: break if i==1: zval=zp else: zval = zval*zp/i i = i + 1 end return ret*2/sqrt(pi) END METHOD
Mean and Variance Calculation
The following method calculates the mean and the biascorrected sample variance of a list of real numbers. It returns a twoitem list containing the mean and that kind of variance in that order. It uses the Welford method presented by J. D. Cook.
METHOD MeanAndVariance(list) if size(list)==0: return [0, 0] if size(list)==1: return [list[0], 0] xm=list[0] xs=0 i=1 while i < len(list) c = list[i] i = i + 1 cxm = (c  xm) xm = xm + cxm *1.0/ i xs = xs + cxm * (c  xm) end return [xm, xs*1.0/(size(list)1)] END