Services & ResourcesWolfram ForumsMathGroup Archive

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

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>