Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Re: Re: Update 6.0.2
  • Next by Date: Re: Pi is not a real number?
  • Previous by thread: Re: Increasing RandomInteger speed
  • Next by thread: Puzzled by the "Variance"