Re: Increasing RandomInteger speed

*To*: mathgroup at smc.vnet.net*Subject*: [mg86568] Re: [mg86490] Increasing RandomInteger speed*From*: Darren Glosemeyer <darreng at wolfram.com>*Date*: Thu, 13 Mar 2008 20:53:33 -0500 (EST)*References*: <200803121026.FAA17301@smc.vnet.net>

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

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