How to generate numbers based on an arbitrary discrete distribution?

Question

How do I generate numbers based on an arbitrary discrete distribution?

For example, I have a set of numbers that I want to generate. Say they are labelled from 1-3 with the follow probabilities:

$$ 1:4\%, 2: 50\%, 3: 46\% $$

The percentages are probabilities that they will appear in the output from the random number generator. I have a pseudorandom number generator that will generate a uniform distribution in the interval [0, 1). Is there any way of doing this?

The number of elements is unbounded, the sum of probabilities is adding up to 1.

Answer 1

One of the best algorithms for sampling from a discrete distribution is the alias method.

The alias method (efficiently) precomputes a two-dimensional data structure to partition a rectangle into areas proportional to the probabilities.

In this schematic from the referenced site, a rectangle of unit height has been partitioned into four kinds of regions--as differentiated by color--in the proportions $1/2$, $1/3$, $1/12$, and $1/12$, in order to sample repeatedly from a discrete distribution with these probabilities. The vertical strips have a constant (unit) width. Each is divided into just one or two pieces. The identities of the pieces and the locations of the vertical divisions are stored in tables accessible via the column index.

The table can be sampled in two simple steps (one for each coordinate) requiring generating just two independent uniform values and $O(1)$ calculation. This improves on the $O(\log(n))$ computation needed to invert the discrete CDF as described in other replies here.

Answer 2

You can do this easily in R, just specify the size you need:

sample(x=c(1,2,3), size=1000, replace=TRUE, prob=c(.04,.50,.46))

Answer 3

In your example, say you draw your pseudorandom Uniform[0,1] value and call it U. Then output:

1 if U<0.04

2 if U>=0.04 and U<0.54

3 if U>=0.54

If the % specified are a, b, ..., simply output

value 1 if U

value 2 if U>=a and U<(a+b)

etc.

Essentially, we are mapping the % into subsets of [0,1], and we know the probability that a uniform random value falls into any range is simply the length of that range. Putting the ranges in order seems the simplest, if not unique, way to do it. This is assuming that you are asking about discrete distributions only; for continuous, can do something like "rejection sampling" (Wikipedia entry).

Answer 4

Suppose there are $m$ possible discrete outcomes. You divide up the interval $[0,1]$ into subintervals based on the cumulative probability mass function, $F$, to give the partitioned $(0,1)$ interval

$$ I_{1} \cup I_{2} \cup \cdots \cup I_{m}$$

where $I_{j} = (F(j-1), F(j))$ and $F(0) \equiv 0$. In your example $m = 3$ and

$$I_1 = (0,.04), \ \ \ \ \ I_2 = (.04,.54), \ \ \ \ \ I_3 = (.54,1)$$

since $F(1) = .04$ and $F(2) = .54$ and $F(3) = 1$.

Then you can generate $X$ with distribution $F$ using the following algorithm:

(1) generate $U \sim {\rm Uniform}(0,1)$

(2) If $U \in I_{j}$, then $X = j$.

• This step can be accomplished by looking at whether $U$ is less than each of the cumulative probabilities, and seeing where the change point (from TRUE to FALSE) occurs, which should be a matter of using a boolean operator in whatever programming language you're using and finding where the first FALSE occurs in the vector.

Note that $U$ will be in exactly one of the intervals $I_{j}$ since they are disjoint and partition $[0,1]$.

Answer 5

One simple algorithm is to start with your uniform random number and in a loop first subtract off the first probability, if the result is negative then you return the first value, if still positive then you go to the next iteration and subtract off the next probability, check if negative, etc.

This is nice in that the number of values/probabilities can be infinite but you only need to calculate the probabilities when you get close to those numbers (for something like generating from a Poisson or negative binomial distribution).

If you have a finite set of probabilities, but will be generating many numbers from them then it could be more efficient to sort the probabilities so that you subtract the largest first, then the 2nd largest next and so forth.

Answer 6

First of all, let me draw your attention to a python library with ready-to-use classes for either integer or floating point random number generation that follow arbitrary distribution.

Generally speaking there are several approaches to this problem. Some are linear in time, but require large memory storage, some run in O(n log(n)) time. Some are optimized for integer numbers and some are defined for circular histograms (for example: generating random time spots during a day). In the above mentioned library I used this paper for integer number cases and this recipe for floating point numbers. It (still) lacks circular histogram support and is generally messy, but it works well.

Answer 7

I had the same problem. Given a set where each item has a probability and whose items' probabilities sum up to one, I wanted to draw a sample efficiently, i.e. without sorting anything and without repeatedly iterating over the set.

The following function draws the lowest of $N$ uniformly distributed random numbers within the interval $[a,1)$. Let $r$ be a random number from $[0,1)$.

\begin{equation} \text{next}(N, a) = 1 - (1 - a) \cdot \sqrt[N]{r} \end{equation}

You can use this function to draw an ascending series $(a_i)$ of $N$ uniformly distributed random numbers in [0,1). Here is an example with $N = 10$:

$a_0 = \text{next}(10, 0)$
$a_1 = \text{next}(9, a_0)$
$a_2 = \text{next}(8, a_1)$
$\dots$
$a_9 = \text{next}(1, a_8)$

While drawing that ascending series $(a_i)$ of uniformly distributed numbers, iterate over the set of probabilities $P$ which represents your arbitraty (yet finite) distribution. Let $0 \leq k < |P|$ be the iterator and $p_k \in P$. After drawing $a_i$, increment $k$ zero or more times until $\sum p_0 \dots p_k > a_i$. Then add $p_k$ to your sample and move on with drawing $a_{i+1}$.

---

Example with the op's set $\{(1, 0.04), (2, 0.5), (3, 0.46)\}$ and sample size $N = 10$:

i  a_i    k  Sum   Draw
0  0.031  0  0.04  1
1  0.200  1  0.54  2
2  0.236  1  0.54  2
3  0.402  1  0.54  2
4  0.488  1  0.54  2
5  0.589  2  1.0   3
6  0.625  2  1.0   3
7  0.638  2  1.0   3
8  0.738  2  1.0   3
9  0.942  2  1.0   3

Sample: $(1, 2, 2, 2, 2, 3, 3, 3, 3, 3)$

---

If you wonder about the $\text{next}$ function: It is the inverse of the probability that one of $N$ uniformly distributed random numbers lies within the interval $[a, x)$ with $x \leq 1$.