Re: Increasing RandomInteger speed

*To*: mathgroup at smc.vnet.net*Subject*: [mg86620] Re: Increasing RandomInteger speed*From*: Mark Fisher <particlefilter at gmail.com>*Date*: Sat, 15 Mar 2008 03:10:34 -0500 (EST)*References*: <200803121026.FAA17301@smc.vnet.net> <frclr0$rob$1@smc.vnet.net>

On Mar 13, 9:53 pm, Darren Glosemeyer <darr... at wolfram.com> wrote: > Art wrote: > > Is there an intelligent way to speed up generating the following > > sample from an inhomogeneous Poisson process to be comparable to the > > homogeneous one? > > > mu = RandomReal[10, 10000]; > > > {t1, sp1} = Timing[RandomInteger[PoissonDistribution[#], 1] & /@ mu]; > > {t2, sp2} = Timing[RandomInteger[PoissonDistribution[10], 10000]]; > > > In[128]:= t1 > > > Out[128]= 8.47253 > > > In[129]:= t2 > > > Out[129]= 0.004 > > > Thanks, > > Art > > The big difference in timing is due to setup time for generation. The > first case has to set up a generator for each of the 10000 elements of > mu, while the second only has to set up a generator once. I can't think > of any way to make the inhomogeneous case as fast as the homogeneous for > small mu values, but it is possible to get values faster by using table > lookup via RandomChoice and weights obtained as a function of the mean. > > This can be done as follows. > > In[1]:= mu = RandomReal[10, 10000]; > > First, determine a cut off for the table. upperbound will cover at least > probability 1 - 10^-10 for each mu<=10. > > In[2]:= upperbound=InverseCDF[PoissonDistribution[10], 1 - 10^-10] > > Out[2]= 36 > > Next, create a rule that maps the probabilities to the values associated > with those probabilities as a fucntion of the mean. In this case, I've > scaled out the Exp[-mm] term because it is a constant for a given mm and > hence we do not need to compute it to get the relative weights (and > computing the exponential for each mean value will add some time to the > set up in the generator we create). > > In[3]:= rl[mm_] = > Table[PDF[PoissonDistribution[mm], x], {x, 0, > upperbound}]/Exp[-mm] -> > Range[0, upperbound]; > > Now Map[]ping RandomChoice with this rule onto mu will give us the > desired sequence of random numbers. > > In[4]:= (vals = Map[RandomChoice[rl[#]] &, mu]); // Timing > > Out[4]= {0.765, Null} > > In[5]:= Timing[sp1 = RandomInteger[PoissonDistribution[#]] & /@ mu;] > > Out[5]= {11.391, Null} > > We can look at moment-based quantities as a rough check of the validity. > > In[6]:= {Mean[#], Variance[#], Skewness[#], Kurtosis[#]} &[vals // N] > > Out[6]= {5.0358, 13.4977, 0.619433, 2.94056} > > In[7]:= {Mean[#], Variance[#], Skewness[#], Kurtosis[#]} &[sp1//N] > > Out[7]= {5.028, 13.2995, 0.628799, 2.99655} > > In cases where the mu values are sufficiently large (with sufficiently > large left as a decision for the user), a normal approximation might be > considered (Poisson approaches a normal for large mu). To do that, > generate n standard normals, scale by the Poisson standard deviations, > shift by the means, and round. This will be much faster because it only > requires normal generation and basic math operations, but it is only an > option when the mean values are all large. > > In[8]:= mu2 = RandomReal[{950, 1050}, 10000]; > > In[9]:= Timing[normapprox = > Round[Sqrt[mu2] RandomReal[NormalDistribution[], Length[mu2]] + > mu2];] > > -16 > Out[9]= {6.37511 10 , Null} > > Darren Glosemeyer > Wolfram Research Following up on Darren's approach, compiling seems to speed things up further: setup[m_, max_] := With[{r = Range[0, max]}, RandomChoice[m^r/r! -> r]] upperBound[n_, eps_] := InverseCDF[PoissonDistribution[n], 1 - eps] cf = Compile[{z}, setup[z, upperBound[10, 10^(-10)]] // Quiet // Evaluate]; mu = RandomReal[10,10^4]; vals1 = cf /@ mu;//Timing vals2 = RandomChoice[rl[#]]& /@ mu;//Timing Through[{Mean, Variance, Skewness, Kurtosis}[# // N]] & /@ {vals1, vals2} // Transpose // Grid // PaddedForm[#, {4, 1}] & --Mark

**References**:**Increasing RandomInteger speed***From:*Art <grenander@gmail.com>