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