# How to generate numbers based on an arbitrary discrete distribution?

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 as follows.

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

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

There are no bounds on how many elements I can have, but the % will add up to 100%.

• I might suggest specifying "...arbitrary discrete distributions" in the title, if that is your question. The continuous case is different. Apr 20, 2012 at 23:56
• A generic way is to perform a binary search within a list of the cumulative probabilities, which in this example would be $(0,0.04,0.54,1.0)$. On the average this takes $\log(n)/2$ probes per generation event. If no probability is extremely small, you can get $O(1)$ performance by creating a vector of equally-spaced values in $[0,1]$ and (in a precomputation stage) assigning an outcome to each value. E.g., in this example you might create the vector $(1,1,1,1,2,\ldots,2,3,\ldots,3)$ (with $50$ 2's and $46$ 3's). Generate a uniform, multiply by 100, and index into this vector: done.
– whuber
Apr 22, 2012 at 21:34
• Also see here Jun 5, 2016 at 22:23
• That "here" link actually links to this very question, @Glen_b... copy-n-paste error? Jun 11, 2016 at 14:50
• @buruzaemon thanks yes that was a mistake; I have corrected it. Jun 12, 2016 at 6:40

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.

• This algorithm is only best if the probabilities are cheap to compute. For example if $n$ is huge it may better not to construct the whole tree. Apr 21, 2012 at 11:19
• +1 So far this is the only reply to suggest and describe an efficient algorithm.
– whuber
May 9, 2012 at 16:01
• Looks a lot like an illustration from Keith Schwarz's page keithschwarz.com/darts-dice-coins If that's the case it would be right to credit that source. Apr 29, 2021 at 0:13

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))

• Personally, i would prefer an algorithm (or somewhere to learn the necessary knowledge), since I am trying to incorporate this into an app that I am building :) Thanks a lot for your answer though :) Apr 20, 2012 at 23:45
• Hmmm ok... Knowing a little more about what you want to do would help us guide you. Can you tell us more about it? (Purpose, context, etc.) Apr 21, 2012 at 0:07
• It is for voting. For example, I have a bunch of photos, and i can only show 6 to a user at a time, i would like to incorporate the "best" to a user at a time, and the user can vote up or down on each photo. The simplest solution that could work right now is the scheme i outlined (each number represent a photo, every down vote would decrease the probability on that photo, and increase on everything else) Apr 21, 2012 at 2:25
• @furtivefelon, you can always port the code from R, o figure out the algorithm from the code and reimplement it. Apr 21, 2012 at 6:49
• I am thinking you might get some good (better) advice on Stack Overflow, since there probably exists some well-known solutions for this specific purpose. I suggest also including the info from your last comment directly into your question. Apr 22, 2012 at 18:15

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

• The algorithm is faster if you sort the categories in decreasing order of probability. That way, you do fewer tests (on average) per random number generated. Apr 21, 2012 at 16:35
• Just to add a quick note on sorting - this will be effective only if you do it once at the start of a sampling scheme - so it won't do well for cases where the probabilities are themselves sampled as part of a larger overall scheme (eg. $p_{j}\sim\text{Dist}$ and then $Pr(Y=j)=p_{j}$). By sorting in this case you are adding the sorting operation into every iteration of sampling - which will be adding $O(n\log(n))$ time to each iteration. However, it may be useful to sort by an approximate guess at the size of the probabilities at the start in this case. Dec 9, 2013 at 11:44

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]$.

• Shouldn't those intervals all be half-closed? Otherwise the boundaries between intervals are not included.. ie. $\{[0,0.04),\ [0.04,0.54),\ [0.54,1]\}$ Apr 21, 2012 at 0:53
• $P(U=u)=0$ for any point $u$ (i.e. the Lebesgue measure of the half open interval is the same as that of the open interval) so I don't think it matters. Apr 21, 2012 at 3:58
• On a finite-precision digital machine, though, maybe someday before the end of the universe it will matter... Apr 21, 2012 at 16:33
• Fair enough, @whuber, see my edit. May 9, 2012 at 16:13
• OK, that is an algorithm. BTW, why don't you just return something like min(which(u < cp))? It would be good to avoid recomputing the cumulative sum on each call, too. With that precomputed, the entire algorithm is reduced to min(which(runif(1) < cp)). Or better, because the OP asks to generate numbers (plural), vectorize it as n<-10; apply(matrix(runif(n),1), 2, function(u) min(which(u < cp))).
– whuber
May 9, 2012 at 16:33

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.

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.

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

• It appears the problem you are addressing abruptly changed in the second paragraph from one that samples from an arbitrary discrete distribution to sampling from a uniform distribution. Its solution appears not to be relevant to the question that was asked here.
– whuber
Jun 11, 2016 at 13:44
• I clarified the last part.
– casi
Jun 13, 2016 at 16:36
• Your answer still seems unrelated to the question. Could you perhaps provide a small but nontrivial worked example of your algorithm? Show us how it would generate a single draw from the set $\{1,2,3\}$ according to the probabilities given in the question.
– whuber
Jun 13, 2016 at 17:03
• I added an example. My answer has something in common with David M Kaplan's answer (stats.stackexchange.com/a/26860/93386), but requires just one instead of N (= sample size) iterations over the set, at the expense of drawing N N-th roots. I profiled both procedures, and mine was much faster.
– casi
Jun 14, 2016 at 16:48
• Thank you for the clarification (+1). It may be of interest to many readers that this isn't a simple random sample, because the outcomes appear in a predetermined, fixed order: a random permutation would have to be applied to the results in order to create a simple random sample. You might also be interested in a parallelizable version of this algorithm in which $$a_j=\frac{\sum_{i=1}^j \log(u_i)}{\sum_{i=1}^{N+1}\log(u_i)}$$ where $u_1,\ldots,u_{N+1}$ is a simple random sample of Uniform(0,1] variates.
– whuber
Jun 14, 2016 at 18:58