Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: easy question about random numbers

• To: mathgroup at smc.vnet.net
• Subject: [mg53516] Re: easy question about random numbers
• From: "Ray Koopman" <koopman at sfu.ca>
• Date: Sat, 15 Jan 2005 01:44:18 -0500 (EST)
• References: <crqb7f$cak$1@smc.vnet.net><200501110631.BAA00684@smc.vnet.net> <cs5bag$3t5$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

DrBob wrote:
> 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: [...]

The alias method, too, runs faster as compiled vector-generating code.
(I apologize for the absence of indentation in what follows, but I'm
posting via Google, which now removes all indentation.)

In[1]:=
<<Statistics`DiscreteDistributions`
<<Graphics`Graphics`

In[3]:=
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]]]

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

Out[7]= {0. Second, Null}
Out[8]= {1.19 Second, -Graphics-}
Out[9]= {0.29 Second, -Graphics-}

In[10]:=
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]}]
Clear@alias
alias[___] := (-1+If[Random[]>c[[#]],a[[#]],#]&)[Random[Integer,{1,n}]]
Clear@alii
alii[samples_] := -1 + Table[(If[Random[] >
c[[#]],a[[#]],#]&)[Random[Integer,{1,n}]],{samples}]
Clear@calii
calii = Compile[{{n,_Integer},{a,_Integer,1},{c,_Real,1},
{samples,_Integer}}, -1 + Table[(If[Random[] >
c[[#]],a[[#]],#]&)[Random[Integer,{1,n}]],{samples}]]

In[17]:=
Timing[{n, a, c} = setal@Table[PDF[dist, i], {i, 0, k}];]
Timing@Histogram@Array[alias, {1*^5}]  (* one at a time *)
Timing@Histogram@alii@1*^5           (*    all at once    *)
Timing@Histogram@calii[n,a,c,1*^5] (* all at once, compiled *)

Out[17]= {0. Second, Null}
Out[18]= {1.74 Second, -Graphics-}
Out[19]= {0.31 Second, -Graphics-}
Out[20]= {0.27 Second, -Graphics-}

• Follow-Ups:
  • Re: Re: easy question about random numbers
    • From: DrBob <drbob@bigfoot.com>

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