Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: easy question about random numbers

• To: mathgroup at smc.vnet.net
• Subject: [mg53474] Re: [mg53424] Re: easy question about random numbers
• From: DrBob <drbob at bigfoot.com>
• Date: Thu, 13 Jan 2005 03:12:25 -0500 (EST)
• References: <crqb7f$cak$1@smc.vnet.net> <200501110631.BAA00684@smc.vnet.net>
• Reply-to: drbob at bigfoot.com
• Sender: owner-wri-mathgroup at wolfram.com

Wow, PieceWise and Compile really do the job; that's quite a bit faster than the alias method, even when ranfun is used to get only one sample at a time. For a Histogram of 100000 samples:

MakeDiscreteRandomFunction[probs_, vals_] := With[{fun = Block[{x},
     Function @@ {x, Piecewise@Transpose[{vals, x
         <= # & /@ (Rest@FoldList[Plus, 0, probs])}]}]},
             Compile[{{n, _Integer}}, fun /@ Table[Random[], {n}]]]
MakeDiscreteRandomFunction[probs_] :=
                           MakeDiscreteRandomFunction[probs,
                             Range[Length[probs]]]

k = 7;
dist = BinomialDistribution[k, 3/10];
Timing[ranfun = MakeDiscreteRandomFunction@Table[PDF[dist, i], {i, 0, k}];]
Timing@Histogram@Flatten@Array[ranfun@1 &, 100000] (* one sample at a time *)
Timing@Histogram@ranfun@100000 (* all the samples at once *)

{0. Second, Null}
{0.5939999999999994*Second, Graphics[]}
{0.1869999999999994*Second, Graphics[]}

setal[p_]:=Block[{n=Length[p],a,b,
     c,i,j,tol},a=Range[n];b=n*p-1;c=Table[1,{n}];tol=2*n*$MachineEpsilon;
     While[{i,j}=Ordering[b][[{1,-1}]];b[[j]]>b[[i]]+
     tol,a[[i]]=j;c[[i]]=1+b[[i]];b[[j]]+=b[[i]];
       b[[i]]=0];{n,a,N[c]}];

Timing[{n, a, c} = setal@Table[PDF[dist, i], {i, 0, k}];]
Clear@alias
alias[___] := (-1 + If[
             Random[] > c[[#1]], a[[#1]], #1] &)[Random[Integer, {1, n}]]
Timing@Histogram@Array[alias, {100000}]

{0. Second, Null}
{0.782*Second, Graphics[]}

Bobby

On Tue, 11 Jan 2005 01:31:12 -0500 (EST), Mark Fisher <mark at markfisher.net> wrote:

> Here's a faster version (again using Piecewise) that incorporates the
> calls to Random[] into the compiled function:
>
> MakeDiscreteRandomFunction[probs_, vals_] :=
>    With[{
>      fun = Block[{x},
>        Function @@ {x, Piecewise @
>          Transpose[{vals,
>            x <= # & /@ (Rest @ FoldList[Plus, 0, probs])}]}
>          ]},
>        Compile[{{n, _Integer}},
>          fun /@ Table[Random[], {n}]
>          ]
>      ]
>
> MakeDiscreteRandomFunction[probs_] :=
>    MakeDiscreteRandomFunction[probs, Range[Length[probs]]]
>
> li2= {1/6, 1/6, 1/6, 1/6, 9/60, 11/60}
> ranfun = MakeDiscreteRandomFunction[li2]
>
> ranfun[10]
>
> --Mark
>
> Pedrito wrote:
>
>> Hi everybody!
>>
>> I wanted to obtain a discrete random number generator that I needed for
>> a project.
>>
>> On the library Statistics`DiscreteDistributions` I could find the DiscreteUniformDistribution
>> function. But I wanted to specify the probability for each one of the
>> states.
>>
>>
>> For instance:
>> If we need to simulate an unfair dice, we could have this probabilities
>> for each one of the sides:
>> {1/6, 1/6, 1/6, 1/6, 9/60, 11/60}
>>
>> So I wrote:
>> li2 = {1/6, 1/6, 1/6, 1/6, 9/60, 11/60}
>> li3=FoldList[Plus,0,li2]
>> Module[{i = 1, r = Random[]}, While[ !li3[[i]] < r < li3[[i + 1]], i++]; i]
>>
>> It works ok but I don't know if there is another (better) way of doing
>> this.
>> Any suggestion?
>>
>
>
>
>



-- 
DrBob at bigfoot.com
www.eclecticdreams.net

• References:
  • Re: easy question about random numbers
    • From: Mark Fisher <mark@markfisher.net>