[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
Prev by Date:
**Re: Problem with transformation rule of a function**
Next by Date:
**Re: Newbie Limit problem**
Previous by thread:
**Re: easy question about random numbers**
Next by thread:
**Re: easy question about random numbers**
| |