In response to a question in comments, here's an outline of a few potentially* faster ways to do discrete distributions than the cdf method.
* I say "potentially" because for some discrete cases a well implemented inverse-cdf approach can be very fast. The general case is harder to make fast without introducing additional tricks.
For the case of four different outcomes as in the example in the question, the naive version of the inverse cdf approach (or effectively equivalent approaches) are fine; but if there are hundreds (or thousands, or millions) of categories it can become slow without being a bit smarter (you certainly don't want to be sequentially searching through the cdf until you find the first category in which the cdf exceeds a random uniform). There are some faster approaches than that.
[You could see the first few things I mention below s having a connection to faster-than-sequential approaches to locating a value using an index and so are in a way just a "smarter version of using the cdf". One can of course look at "standard" approaches to solve related problems like "searching a sorted file", and end up with methods with much faster than sequential performance; if you can call suitable functions, such standard approaches may often be all you need.]
Anyway, to some efficient approaches for generating from discrete distributions.
1) the "Table method". Instead of being $O(n)$ for $n$ categories, once set up, the "simple" version of this in (a) (if the distribution is suitable) is $O(1)$.
a) Simple approach - assuming rational probabilities (done on the above data example):
- set up an array with 10 cells, containing a '1', four '2's, two '3's and three '4's. Sample that using a discrete uniform (easy to do from a continuous uniform), and you get simple fast code.
b) More complex case - doesn't need 'nice' probabilities. Use $2^k$ cells, or rather, you'll end up using a few less than that. So for example, consider the following:
x 0 1 2 3 4 5 6
P(X=x) 0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002
(We could have 10000 cells and use the previous exact approach, of course, but what if these probabilities are irrational, say?)
Let's use $k=8$. Multiply the probabilities by $2^k$ and truncate to find out how many of each cell type we need:
x 0 1 2 3 4 5 6 Tot
P(X=x) 0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002 1.0000
[256p(x)] 117 0 50 84 0 2 0 253
Then the last 3 cells are basically "generate instead from this other distribution" (i.e. p(x) - \frac{\lfloor 256 p(x)\rfloor}{256} normalized to a p.m.f.):
x* 0 1 2 3 4 5 6 Tot
P(X=x*) 0.091200 0.273067 0.272000 0.142933 0.187733 0.016000 0.017067 1.000000
The "spillover" table can be done by any reasonable method (you only get here about 1% of the time, it doesn't need to be as fast). So $\frac{253}{256}$ of the time we generate a random uniform, use its first 8 bits to pick a random cell, and output the value in the cell; after the initial setup all of this can be made very fast. The other $\frac{3}{256}$ of the time we hit a cell that says "generate from the second table". Almost always, you generate a
single uniform on $(0,1)$ and get a discrete random number from a multiply, a truncation and the cost of accessing an array element.
2) "Squaring the histogram" method; this is kind of related to (1), but each cell can actually generate either one of two values, depending on a (continuous) uniform. So you generate a discrete value from 1 to n, then within each one, check whether to generate its main value or its second value. It works with bounded random variables. There's no spillover table, and it generally uses much smaller tables than method (1). Usually, it's set up so that the choice of 1:n uses the first few bits of a uniform random number, and the remainder of it tells you which of the two values for that bin to output.
Perhaps the easiest way to outline the method is to do it on the above example:
Think of the distribution as a histogram with 4 bins:
We cut off the tops of the tallest bars and put them in the shorter ones, 'squaring it off'. The average 'height' of a bar will be 0.25. So we cut 0.15 off the second bar and put it in the first and 0.05 off the fourth and put it in the third:
Its always possible to organize this in such a way that no bin ends up with more than 2 colors, though one color may end up in several bins.
So now you choose one of the 4 bins at random (requires 2 random bits off the top of a uniform). You then use the remaining bits to specify a uniformly distributed vertical position and compare with the break between colors to work out which of two values to output. While very fast it's usually not quite as fast as the 'table' method.
--
These methods can be adapted to deal with unbounded variables, where again, it's 'mostly fast'.
A reference: http://www.jstatsoft.org/v11/i03/paper
The relatively slow part of these is creating the tables of values; they're suitable when you know what you're going to be generating ("we need to sample values from this distribution many times in the future") rather than trying to create it as you go. "We need to sample a million values from this ASAP, but we'll never need to do it again" creates different priorities; in many situations some of the "standard computing approaches" to looking up sorted values (i.e. to doing the cdf method more quickly) may actually be about your best choice.
There are still other fast approaches to generating from discrete distributions. Carefully coded, you can do some very fast generation. For example:
3) the rejection method ("accept-reject") can be done with discrete distributions; if you have a discrete majorizing function ("envelope") which is a scaled-up discrete pmf that you can already generate from in a fast way, it adapts directly, and in some cases can be very fast. More generally you can take advantage of being able to generate from continuous distributions (for example by discretizing the result back to a discrete envelope).
Here imagine that we have some discrete probability function $f$ for which we don't have a convenient cdf (or inverse-cdf) -- indeed in this illustration we didn't even have the normalizing constant, so our plot is unnormalized:
Now we need to find some convenient-to-generate-from discrete probability function $g$, which can be multiplied by a constant $c$ and be everywhere at least as large as $f$ (we need to be sure that this remains true for all $x$ values). That is, $c.g(x)\geq f(x)$ for all the possible $x$ values.
Sometimes a suitable $g$ is easily identified, but one useful option is to take a mixture of a discrete uniform for the left part, and a distribution with at least as heavy a tail as $f$ on the right. Two reasonably convenient choices for that are a geometric distribution (when the tail isn't decreasing slower than exponentially) and something like a discretized Pareto or discretized half-Cauchy distribution, obtained by taking $\lfloor X\rfloor$ for some Pareto or half-Cauchy random variate $X$ (in either case for when the pmf is decreasing slower than exponentially).
(For that matter, the geometric itself can be generated by discretizing an exponential.)
In this case, a discrete uniform on the left and a geometric on the right work quite well:
(Reminder: what is plotted here is an unnormalized pmf, so the y-axis doesn't represent probability but something proportional to probability)
Then the procedure is a matter of simulating a proposed value $x$ from $g$, simulating a uniform, $U$ on $(0,c.g(x))$ and if $U<f$, accepting that proposed $x$ (otherwise rejecting it and generating a new proposed $x$).
