Re: Re: newbie is looking for a customDistribution function

*To*: mathgroup at smc.vnet.net*Subject*: [mg50460] Re: [mg50435] Re: newbie is looking for a customDistribution function*From*: DrBob <drbob at bigfoot.com>*Date*: Sat, 4 Sep 2004 01:44:01 -0400 (EDT)*References*: <200409032122.i83LM3O2023123@rm-rstar.sfu.ca>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

I'm sure the Alias Method is brilliant, but I'd like to see an implementation. Here's a decent method---NOT for the OP's original problem, but for a more general case where discrete probabilities are known. n is the number of discrete values and p is the vector of their probabilities. "inverse" is the inverse CDF. n = 300; p = (#1/Tr[#1] &)[(Random[] &) /@ Range[n]]; inverse = Interpolation[Transpose@{FoldList[Plus, 0, p], Range[0, Length@p]}, InterpolationOrder -> 0]; random2[] := inverse@Random[] samples = 1000000; result = Replace[MapAt[Frequencies, Timing@Array[random2[] &, samples], 2], {a_, b_} :> {N[ a/samples], b}, 2]; First@result 4.343 Second result[[-1, All, 1]] - p // Abs // Max 0.000186455 The last number is the maximum absolute difference between observed frequencies and the corresponding p values. (That statement fails if not all discrete values are represented in the sample.) Timings on my machine for a million samples at n=3, 30, and 300 (the number of discrete values) are between 4.3 and 4.4 seconds for all three values of n. Bobby On Fri, 03 Sep 2004 14:22:03 -0700, Ray Koopman <koopman at sfu.ca> wrote: > On Fri, 03 Sep 2004 10:15:56 -0500 drbob at bigfoot.com wrote: > >> On Fri, 3 Sep 2004 03:36:14 -0400 (EDT), Ray Koopman <koopman at sfu.ca> > >> wrote: > >>> However, if you really want to use the distribution instead of the data > >>> that gave rise to it then you should look into the "Alias Method" of > >>> generating random observations from an arbitary discrete distribution. > >> > >> How would we look into that? Is that in the Mathematica help files, > >> math books, what? > > > > Before I wrote that, I Googled for <"alias method" "random number"> > > and got ~191 references. I assumed everyone would do something similar. > > In retrospect, I should have just pointed to > > > > http://cgm.cs.mcgill.ca/~luc/rnbookindex.html > > > > which is a source that everyone who uses anything more complicated than > > simple Uniform random numbers should know about. > > > > Ray > > -- DrBob at bigfoot.com www.eclecticdreams.net