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.

  • 1
    $\begingroup$ RandomVariate[MultivariatePoissonDistribution[1, m-1]] is equally slow, in case anyone was about to try. $\endgroup$
    – ssch
    Nov 5, 2013 at 17:32

4 Answers 4


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)

  • $\begingroup$ Thanks: Compiling to C speeds it up a bit more. $\endgroup$
    – asim
    Nov 5, 2013 at 22:33
  • $\begingroup$ @asim This is exactly what I wanted to implement! However, I posted my own version in update in my answer. $\endgroup$
    – ybeltukov
    Nov 5, 2013 at 23:38
  • 3
    $\begingroup$ Related: imgs.xkcd.com/comics/poisson.jpg :) $\endgroup$
    – ybeltukov
    Nov 6, 2013 at 1:04

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.

  • 1
    $\begingroup$ +1. Perhaps it's more high level to do EmpiricalDistribution[PDF[ ParameterMixtureDistribution[PoissonDistribution[par], par \[Distributed] UniformDistribution[{1, 5}]], zs] -> zs] $\endgroup$
    – Rojo
    Nov 5, 2013 at 18:18
  • $\begingroup$ @Rojo I was looking for a such function with no luck. Thanks! $\endgroup$
    – ybeltukov
    Nov 5, 2013 at 18:24
  • 1
    $\begingroup$ @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. $\endgroup$
    – asim
    Nov 5, 2013 at 19:00
  • 1
    $\begingroup$ @asim Now I understand you! I will think how to deal with it. $\endgroup$
    – ybeltukov
    Nov 5, 2013 at 19:09
  • 1
    $\begingroup$ A 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. $\endgroup$ Mar 15, 2019 at 13:44

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.

  • $\begingroup$ 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. $\endgroup$
    – asim
    Nov 5, 2013 at 17:11
  • $\begingroup$ @asim Is the Julia generator an exact generator ... or does it use Gaussian or other approximation? $\endgroup$
    – wolfies
    Nov 5, 2013 at 17:29
  • $\begingroup$ I think it relies on the Rmath standalone c library that underlies the distribution functionality of the R statistical language. $\endgroup$
    – asim
    Nov 5, 2013 at 17:32

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]]
  • $\begingroup$ +1 My tests show that it is fastest solution! (with options RuntimeAttributes -> {Listable}, Parallelization -> True, CompilationTarget -> "C") $\endgroup$
    – ybeltukov
    Nov 6, 2013 at 12:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.