> osama wrote: > > Hallo; > Do you know of a fast algorithm to generate random numbers from a > discrete distribution? I am looking for a disrtibution with pmf =x^a > p^x/c
Presumably for x = 1, 2, 3, ....? Or on some other set of x values? As a first rough attempt, I'd probably use rejection based on rounding or even just rounding up a gamma r.v. (note how a relates to (alpha-1) and p to exp(-1/beta) in the gamma - those or very similar values are probably what you want). For example, imagine your pmf is replaced by a continuous density of piecewise constant height (i.e. with height x*^a p^x*/c, where x* = round(x) - that is, by a step function/"histogram" -type of density. Note that since the width of wach bar is 1, this still has total probability 1. Now draw over it a multiple of the gamma density that has alpha = a+1 and beta = -1/ln(p), with the multiple chosen so that the height of the curve is always at least as high as the step density above (a little analysis required to find where the ratio is at a minimum, which determines your constant (another possibility is to use a different density to the left and right of the mode). Now do standard accept/reject, with low probability of rejection. In practice, I'd probably set up one of George Marsaglia's quick methods for generating from discrete distributions, with the extreme right tail generated by accept/reject. (Table lookup is the method I'd most likely use, but squared histogram is also quite good). For a brief description of both of those, see the thread "Simulating from a Histogram" in sci.math.num-analysis (posts between 5 and 18 March 1999). CHeck it out on groups.google.com. I found it by searching for: george marsaglia fast discrete random and it came up as the first item. George's contributions are the 5th and 9th posts to the thread and bear careful reading. George's description of the square histogram method in his second is far more illuminating than Knuth's discussion of the related "alias" method. [Note that if you use table lookup, there's *two* leftover parts you need to generate - the tops of each bit of pmf not covered by the table, as well as the extreme right tail.] Glen . . ================================================================= Instructions for joining and leaving this list, remarks about the problem of INAPPROPRIATE MESSAGES, and archives are available at: . http://jse.stat.ncsu.edu/ . =================================================================
