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>
• 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

```

• Prev by Date: Re: RE: ExpandAll Problem with Rules
• Next by Date: Re: Parallel Toolkit Example
• Previous by thread: Re: Re: newbie is looking for a customDistribution function
• Next by thread: Re: newbie is looking for a customDistribution function