Why is Poisson Random Deviate Generation so slow?

Question

I am generating Poisson deviates for some numerical work. Mathematica 9.0.1 is very slow in generating these random numbers, as can be seen below.

In[1]:= nobs = 100000;

In[2]:= AbsoluteTiming[
 m = RandomReal[{1, 5}, nobs];
 m1 = Map[RandomVariate[PoissonDistribution[#]] &, m];
 ]

Out[2]= {10.194500, Null}

In[3]:= Mean[m1] // N

Out[3]= 2.99884

I need fast Poisson deviates for any vector m (i.e., it does not need to contain Uniform draws). The above is very slow for serious numerical work. A comparable function in Julia, is significantly faster, as can be seen below.

function pois(nobs::Integer)
         m1=rand(Uniform(1,5), nobs)
         for i in 1:nobs
            m1[i]=rand(Poisson(m1[i]))
         end
         m1
end

In [16]: @time draws=pois(nobs);

elapsed time: 0.005883954 seconds (800200 bytes allocated)

Is there any way one can speed up the Mathematica Poisson random number generation? Compilation does not work, unfortunately, and none of the distribution functionality is vectorized automatically.

Hopefully, future versions will address this aspect.

Answer 1

Here is a standard algorithm to generate poisson-distributed numbers that is probably fast enough, but it is not advised to use this for large values of λ as the exponential is prone to inaccuracy, so for anything greater than 7,8,9,10 (seems to be no consensus on "large", perhaps because we have larger and larger floats) you're probably better of to use normal-distribution approximation.

simplePoisson = Compile[{{λ, _Real}},
  Module[{b = 1., i, a = Exp[-λ]},
   For[i = 0, b >= a, i++, b *= RandomReal[]];
   i - 1
   ],
  RuntimeAttributes -> {Listable},
  Parallelization -> True]

nobs = 100000;
m = RandomReal[{1, 5}, nobs];
m1 = simplePoisson[m]; // AbsoluteTiming
(* 0.04s *)

Sanity check, generate a bunch of samples from $Poi(5)$:

m1 = simplePoisson[ConstantArray[5, 10000]];
GraphicsRow[{
  ProbabilityPlot[m1, PoissonDistribution[5]],
  QuantilePlot[m1, PoissonDistribution[5], PlotRange -> All]}]

It looks good thankfully.

Computer generation of Statistical Distributions (pdf)

Answer 2

To obtain the resulting distribution you need to integrate the Poisson distribution with the uniform distribution

a = 1;
b = 5;
p = Integrate[(E^-y y^z)/z!, {y, a, b}]/(b - a)

(Gamma[1 + z, 1] - Gamma[1 + z, 5])/(4 z!)

Now the easiest way is to construct empirical distribution from numerical values of probabilities

maxz = 40;
d = EmpiricalDistribution[Table[N[p], {z, 0, maxz}] -> Range[0, maxz]];

Probability of 40 is less then the machine precision so maxz = 40 is enough.

AbsoluteTiming[res = RandomVariate[d, 100000];]

{0.008653, Null}

It is really fast!

Verification:

Histogram[{res, m1}, {0, 15, 0.5}]

Also you can calculate Mean, Variance, etc. directly from the distribution

Mean[d]

3.

Variance[d]

4.33333

Update

For fixed set of $\lambda$ you can use the following one-liner

λ = RandomReal[{1, 5}, 100000];

res = Total[#,{2}]&@UnitStep[λ + Transpose@Accumulate@Log@RandomReal[1, {14, Length[λ]}]];

Here 14 is the maximum expected number. Timing is ~2 times slower then ssch's compiled solution. I think it is not bad for such simple code!

To understand this code just take the logarithm of ssch's code.

Answer 3

Mathematica's in-built generators work best when you have a common distribution, say Poisson(4), and you want to generate 100,000 pseudo-random drawings from it. But, your example does not have a common parent ... in fact, every one of your drawings is from a different parent: 1 drawing from Poisson(2.2), 1 drawing from Poisson(4.3), 1 drawing from Poisson(3.7) etc ... which is highly inefficient. I don't know the quality of the Julia generator you refer to, or if it matches the precision of the Mma generator.

In any event, you can easily speed up your code by simply parallelising it using ParallelMap instead of Map:

ParallelMap[RandomVariate[PoissonDistribution[#]] &, m]; // AbsoluteTiming

{2.497513, Null}

which is about 6 times faster on my Mac than your original:

Map[RandomVariate[PoissonDistribution[#]] &, m]; // AbsoluteTiming

{16.601512, Null}

---

As an aside, ... trying to manually 'parallelise' random number generation in Mma usually just slows things down: this is because (i) the in-built functions appear to take advantage of multiple cores automatically anyway, and (ii) forcing manual parallel code just ends up shifting huge amounts of data between the various kernels. BUT, here, in your example, it actually speeds things up beautifully ... precisely because every one of your drawings is from a different Poisson parent.

Answer 4

I've been using one form or another of the following Poisson generator since v3.
It needs only one random value but does a little more arithmetic inside the loop.

genpoi = Compile[{{m,_Real}},
Module[{t = 1.-RandomReal[]*E^m, x = 0, mi = m^-1},
While[Negative@t, t = ++x*t*mi + 1.]; x]]