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>