Re: Re: Increasing RandomInteger speed
- To: mathgroup at smc.vnet.net
- Subject: [mg86615] Re: [mg86568] Re: [mg86490] Increasing RandomInteger speed
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Sat, 15 Mar 2008 03:09:37 -0500 (EST)
- References: <200803121026.FAA17301@smc.vnet.net> <18813851.1205537280821.JavaMail.root@m08>
- Reply-to: drmajorbob at longhorns.com
Just for grins, here's a bit more comparison of the two samples.
First to duplicate Darren's work:
mu = RandomReal[10, 10000];
upperbound = InverseCDF[PoissonDistribution[10], 1 - 10^-10]
rl[mm_] =
Table[PDF[PoissonDistribution[mm], x], {x, 0, upperbound}]/
Exp[-mm] -> Range[0, upperbound];
(vals = Map[RandomChoice[rl[#]] &, mu]); // Timing
Timing[sp1 = RandomInteger[PoissonDistribution[#]] & /@ mu;]
{Mean[#], Variance[#], Skewness[#], Kurtosis[#]} &[vals // N]
{Mean[#], Variance[#], Skewness[#], Kurtosis[#]} &[sp1 // N]
36
{0.394639, Null}
{5.99402, Null}
{4.9936, 13.4977, 0.629637, 2.93044}
{4.9787, 13.5706, 0.636508, 2.9522}
I think these stats can be close for significantly different CDF's, so.. .
Clear[pairs, cdf]
pairs[x_List] :=
Last /@ Split[Transpose@{Sort@x, Range@Length@x/Length@x},
First@#1 == First@#2 &]
cdf[x_List] :=
cdf[x] = Interpolation[pairs@x, InterpolationOrder -> 1]
one = Show[ListPlot@pairs@vals, Plot[cdf[vals]@k, {k, 0, upperbound}],
PlotRange -> All];
two = Show[ListPlot@pairs@sp1, Plot[cdf[sp1]@k, {k, 0, upperbound}],
PlotRange -> All];
Show[one, two]
Plot[cdf[vals]@k - cdf[sp1]@k, {k, 0, 36}]
BTW, I really like the fact that interpolation functions default to
error-free extrapolation now. I'm not sure when that was changed, but I
much prefer it.
Bobby
On Thu, 13 Mar 2008 20:53:33 -0500, Darren Glosemeyer
<darreng 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
>
>
--
DrMajorBob at longhorns.com
- References:
- Increasing RandomInteger speed
- From: Art <grenander@gmail.com>
- Increasing RandomInteger speed