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