Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • 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