Why is Poisson Random Deviate Generation so slow?



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

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.


Posted 2013-11-05T16:57:24.587

Reputation: 1 855

1RandomVariate[MultivariatePoissonDistribution[1, m-1]] is equally slow, in case anyone was about to try. – ssch – 2013-11-05T17:32:27.743



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]];
  ProbabilityPlot[m1, PoissonDistribution[5]],
  QuantilePlot[m1, PoissonDistribution[5], PlotRange -> All]}]


It looks good thankfully.

Computer generation of Statistical Distributions (pdf)


Posted 2013-11-05T16:57:24.587

Reputation: 16 150

Thanks: Compiling to C speeds it up a bit more. – asim – 2013-11-05T22:33:24.930

@asim This is exactly what I wanted to implement! However, I posted my own version in update in my answer. – ybeltukov – 2013-11-05T23:38:53.563


Related: http://imgs.xkcd.com/comics/poisson.jpg :)

– ybeltukov – 2013-11-06T01:04:46.600


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!


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

enter image description here

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






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.


Posted 2013-11-05T16:57:24.587

Reputation: 41 907

1+1. Perhaps it's more high level to do EmpiricalDistribution[PDF[ ParameterMixtureDistribution[PoissonDistribution[par], par \[Distributed] UniformDistribution[{1, 5}]], zs] -> zs] – Rojo – 2013-11-05T18:18:27.493

@Rojo I was looking for a such function with no luck. Thanks! – ybeltukov – 2013-11-05T18:24:03.950

@ybeltukov This is good, but I used the uniform just to have different Poisson rate parameters on each draw of the Poisson. I am not interested in a mixture of a poisson and the uniform. – asim – 2013-11-05T18:49:29.093

@asim It is not a mixture distribution. I think it is exactly what you need. My res and your m1 has identical statistical properties.

– ybeltukov – 2013-11-05T18:56:32.973

1@ybeltukov You did a great job, but I need a fast Poisson for any m vector that includes positive numbers. I am not interested, particularly, in using those that are generated from a Uniform distribution. – asim – 2013-11-05T19:00:42.760

1@asim Now I understand you! I will think how to deal with it. – ybeltukov – 2013-11-05T19:09:02.260

@ybeltukov Your first construction is a parameter-mix distribution ... and very nice too ... though I am not sure that deriving the symbolic parameter-mix distribution retains the desired intention to generate Monte Carlo drawings of BLAH, because part of that process has now become theoretical/exact rather than simulated. – wolfies – 2013-11-06T13:03:07.057

1A more compact and more stable way to express p is GammaRegularized[1 + z, a, b]/(b - a); this avoids subtractive cancellation for large values of z. – J. M.'s ennui – 2019-03-15T13:44:38.203


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.


Posted 2013-11-05T16:57:24.587

Reputation: 8 048

I am aware that parallelization can speed up things. But still, even after this, it is too slow, in the context of numerical work involving thousands of iterations. It would be great if the function could be compiled, at least. – asim – 2013-11-05T17:11:45.297

@asim Is the Julia generator an exact generator ... or does it use Gaussian or other approximation? – wolfies – 2013-11-05T17:29:02.140

I think it relies on the Rmath standalone c library that underlies the distribution functionality of the R statistical language. – asim – 2013-11-05T17:32:04.897


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]]

Ray Koopman

Posted 2013-11-05T16:57:24.587

Reputation: 3 196

+1 My tests show that it is fastest solution! (with options RuntimeAttributes -> {Listable}, Parallelization -> True, CompilationTarget -> "C") – ybeltukov – 2013-11-06T12:06:16.050