Re: easy question about random numbers

*To*: mathgroup at smc.vnet.net*Subject*: [mg53400] Re: [mg53382] easy question about random numbers*From*: DrBob <drbob at bigfoot.com>*Date*: Sun, 9 Jan 2005 23:04:06 -0500 (EST)*References*: <200501090402.XAA12446@smc.vnet.net> <opskbs88l9iz9bcq@monster.ma.dl.cox.net> <courier.41E1CD8F.00006418@softhome.net>*Sender*: owner-wri-mathgroup at wolfram.com

I don't understand the alias method either, frankly... but it works! Your method is fast and simple when it works, but here's an easy case where it doesn't: li2 = (#1/Plus @@ #1 & )[Table[Log[Random[]], {6}]] li4 = Rationalize[li2*LCM @@ Denominator[Rationalize[li2]]] a = Flatten[Table[Table[i, {li4[[i]]}], {i, Length[li4]}]] {0.15114701195509367, 0.03837478964581514, 0.056101782183962834, 0.4218879106253158, 0.07794091287875893, 0.25454759271105365} {0.15114701195509367, 0.03837478964581514, 0.056101782183962834, 0.4218879106253158, 0.07794091287875893, 0.25454759271105365} {} Maybe you were intending to handle only rational probabilities, but here's an example with rational probabilities that also fails (not enough memory, so the kernel shuts down): n = 8; dist = BinomialDistribution[n, 3/10]; li2 = Table[PDF[dist, i], {i, 0, n}]; li4 = Rationalize[li2 LCM @@ Denominator[Rationalize[li2]]]; a = Flatten@Table[Table[i, {li4[[i]]}], {i, Length[li4]}]; If I replace 8 with 7 in that code, look how long a turns out to be: n=7; dist=BinomialDistribution[n,3/10]; li2=Table[PDF[dist,i],{i,0,n}]; li4=Rationalize[li2 LCM@@Denominator[Rationalize[li2]]]; a=Flatten@Table[Table[i,{li4[[i]]}],{i,Length[li4]}]; Length@a 10000000 The alias method has no problem with these examples. For instance: n = 8; dist = BinomialDistribution[n, 3/10]; li2 = Table[PDF[dist, i], {i, 0, n}]; Timing[{n, a, c} = setal[li2]] Clear@alias alias[___] := (-1 + If[Random[] > c[[#1]], a[[#1]], #1] &)[Random[Integer, {1,n}]] Timing@Histogram@Array[alias, {100000}] {0. Second, {9, {4, 3, 5, 5, 5, 3, 2, 4, 3}, {0.518832, 0.86887, 0.957819, 0.816951, 1., 0.420079, 0.0900169, 0.0110225, 0.00059049}}} {0.7350000000000002*Second, Graphics[]} Bobby On Sun, 09 Jan 2005 17:34:23 -0700, Pedrito <pedrito6 at softhome.net> wrote: > Thanks for helping me, Bobby. > > I have been trying the differents methods for creating random numbers and > I think the faster are these two: > > -------------- First method > li2 = {1/6, 1/6, 1/6, 1/6, 9/60, 11/60}; > li4 = Rationalize[li2 LCM @@ Denominator[Rationalize[li2]]]; > a = Flatten[Table[Table[i, {li4[[i]]}], {i, Length[li4]}]]; > > And now it's necessary to obtain by hand the length of "a": > Length[a] > Out[]:=60 > > Finally, this is the function: > (a[[Random[Integer, {1, 60}]]] & )[] > > -------------- Second method > 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]}]; > p = (#1/Plus @@ #1 & )[Table[Log[Random[]], {6}]]; {n, a, c} = setal[p]; > > And this is the function: > (If[Random[] > c[[#1]], a[[#1]], #1] & )[Random[Integer, {1, n}]] > > -------------------------------------- > > But I'm sorry to say that I'm not able of understanding the second > method (the "alias" method). > What do you think about the first method? > > > Pedrito > > > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**References**:**easy question about random numbers***From:*Pedrito <pedrito6@softhome.net>