Seeking reference for O(1) random variable lookup algorithm
me13013
free to suggest more appropriate newsgroups.
There is a clever algorithm stuck in my memory, for which I have
forgotten the reference (and the author to credit). I have had no
success searching for a reference online either. I am hoping someone
here can point me to a reference. I know a reference exists, because
I remember discovering the algorithm in some journal paper, probably
written in the 60s or 70s. I didn't stumble upon it until about 2000.
The algorithm solves the problem of picking a random variable from a
finite number of classes, with the probability of each class being
arbitrary (but given, and constant over time. It is able to do a
lookup in constant time, after some pre-computation to create a lookup
table.
How it works is best described by an example. Suppose I have ten
classes and the probabilities are these:
p(A) = .077
p(B) = .093
p(C) = .044
p(D) = .091
p(E) = .126
p(F) = .147
p(G) = .016
p(H) = .168
p(I) = .169
p(J) = .069
It builds a table with 10 entries (shown below). Conceptually each
entry is a unit bar, divided into two pieces (the dividing points are
generally different for each bar). The left piece is assigned to one
class, the right to another (some bars may have both pieces assigned
to the same class). To convert a unit random value r to a random
selection, we multiply r by the number of classes, lookup the
corresponding table entry (using as index the floor of the scaled r).
Then we compare the scaled r to the entry's split value to decide
which side of the bar we're on, and which symbol to emit.
For example, r=.503 scales to 5.03; entry [5] is "D 5.91 H", and
since 5.03 is less than 5.91, we select D.
[0] G 0.16 I
[1] C 1.44 H
[2] J 2.69 F
[3] A 3.77 E
[4] I 4.85 F
[5] D 5.91 H
[6] B 6.93 H
[7] H 7.96 E
[8] E 8.99 F
[9] F 9.50 F
Table construction is not difficult but I'll leave that unexplained
unless anyone wants to see it.
Does this algorithm ring a bell?
Thanks for any help,
Bob H
Ilmari Karonen
>
> It builds a table with 10 entries (shown below). Conceptually each
> entry is a unit bar, divided into two pieces (the dividing points are
> generally different for each bar). The left piece is assigned to one
> class, the right to another (some bars may have both pieces assigned
> to the same class). To convert a unit random value r to a random
> selection, we multiply r by the number of classes, lookup the
> corresponding table entry (using as index the floor of the scaled r).
> Then we compare the scaled r to the entry's split value to decide
> which side of the bar we're on, and which symbol to emit.
>
> unless anyone wants to see it.
>
> Does this algorithm ring a bell?
Nope, sorry, but it's a clever algorithm and I'll be sure to keep it
mind. Thanks!
I can see how to construct the table in O(n log n) time using a
min-max heap. (Pick the classes with the lowest and the highest
probability, remove the former, reduce the probability of the latter
by 1/n minus the probability of the former, and repeat.) Is there a
faster and/or a simpler way?
Note that it's easy to force the table size to be, say, a power of 2
simply by introducing some dummy classes with zero probabilities.
--
Ilmari Karonen
To reply by e-mail, please replace ".invalid" with ".net" in address.
me13013
> Nope, sorry, but it's a clever algorithm and I'll be sure to keep it
> mind. Thanks!
I haven't been able to find this anywhere on the web. I'm surprised
it is not better known.
I have a faint idea that this algorithm is due to Knuth (but I could
be wrong). Last night I was trying to retrace the steps that led me
to it ten years ago. I think I was looking at existing methods to
generate RV from a normal distribution, and found a survey paper that
mentioned other sampling methods, which pointed me at this one.
Attempts to recreate that search using google scholar were not
fruitful.
> I can see how to construct the table in O(n log n) time using a
> min-max heap. (Pick the classes with the lowest and the highest
> probability, remove the former, reduce the probability of the latter
> by 1/n minus the probability of the former, and repeat.) Is there a
> faster and/or a simpler way?
What you describe is (almost) exactly how I do it. I deal with the
probabilities scaled up by N, so your 1/n becomes a simpler 1.
I'm not sure if there is a faster way. Of course if the number of
classes is moderate and the number of samples to be generated is
large, the pre-computation time probably doesn't matter.
> Note that it's easy to force the table size to be, say, a power of 2
> simply by introducing some dummy classes with zero probabilities.
The advantage being to speed up the multiplication that produces the
table index. Nice, I hadn't thought of that.
Another minor improvement, if the classes are just the integers
0..N-1, is the order the table rows by the left symbol and omit the
left symbols from table.
Bob H
Ilmari Karonen
It also simplifies things if you're working directly with integers.
For example, if your RNG gives you n-bit unsigned integers, you can
just use the k highest bits as an index to the table and compare the
n-k lowest bits against the (pre-scaled) threshold to decide which
class to pick.
(Or you could do it the other way around and save one shift, but I
wouldn't recommend trying that with a LCRNG, since their low bits are
not as random as the high bits. With better generators it should be
OK, though, but a shift is pretty cheap in any case.)
> Another minor improvement, if the classes are just the integers
> 0..N-1, is the order the table rows by the left symbol and omit the
> left symbols from table.
Hmm, true. And you don't even need an extra sorting pass, since you
can construct the table that way to begin with. Clever.
me13013
> On 2010-01-14, me13013 <me13...@gmail.com> wrote:
> >> On Jan 14, 5:54 am, Ilmari Karonen <usen...@vyznev.invalid> wrote:
> >> Note that it's easy to force the table size to be, say, a power of 2
> >> simply by introducing some dummy classes with zero probabilities.
>
> > The advantage being to speed up the multiplication that produces the
> > table index. Nice, I hadn't thought of that.
>
> It also simplifies things if you're working directly with integers.
> For example, if your RNG gives you n-bit unsigned integers, you can
Yep, exactly. In the past, the way I implemented the index operation
was, given 32-bit r, multiply by n in a 64-bit register, then use the
most significant 32 bits as the index. Adding the extra table entries
allows me to replace the multiply with a right shift.
Bob H
geo
See "Fast generation of discrete random variables",
G. Marsaglia, W.W Tsang and J. Wang
Journal of Statistical Software, July 2004,v11, Issue 3
based of the original of Marsaglia:
"Generating Discrete Random Variables in a Computer",
Comm. ACM, 6, 37-38, 1963
me13013
> See "Fast generation of discrete random variables",
> G. Marsaglia, W.W Tsang and J. Wang
> Journal of Statistical Software, July 2004,v11, Issue 3
> based of the original of Marsaglia:
> "Generating Discrete Random Variables in a Computer",
> Comm. ACM, 6, 37-38, 1963
Thanks, George. Looking at those led me to what I think is the
reference that I was thinking of. Which is
A. J. Walker, “An efficient method for generating discrete
random variables with general distributions,” ACM
Transactions on Mathematical Software, vol. 3, pp.
253–256, 1977.
It is interesting that I took something of the same path I took back
around 2000. I first found your 1998 paper on the Monty Python
method, and remembered the diagrams therein. That led me to the 1986
DeVroye treatise, which in turn led me to Walker.
Your 2004 paper improves on this method by adding a very fast front
end. Definitely worth reading.
Thanks,
Bob H