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

MathGroup Archive 2008

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

Search the Archive

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


  • Prev by Date: Re: Unexpected failures of FullSimplify on some Root elements
  • Next by Date: Re: Re: Re: Update 6.0.2
  • Previous by thread: Increasing RandomInteger speed
  • Next by thread: Re: Increasing RandomInteger speed