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>
- Increasing RandomInteger speed