# 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:= nobs = 100000;

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

Out= {10.194500, Null}

In:= Mean[m1] // N

Out= 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 : @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.

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

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],
QuantilePlot[m1, PoissonDistribution, PlotRange -> All]}] It looks good thankfully.

Computer generation of Statistical Distributions (pdf)

• Thanks: Compiling to C speeds it up a bit more. – asim Nov 5 '13 at 22:33
• @asim This is exactly what I wanted to implement! However, I posted my own version in update in my answer. – ybeltukov Nov 5 '13 at 23:38
• Related: imgs.xkcd.com/comics/poisson.jpg :) – ybeltukov Nov 6 '13 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!

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.

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

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

• +1 My tests show that it is fastest solution! (with options RuntimeAttributes -> {Listable}, Parallelization -> True, CompilationTarget -> "C") – ybeltukov Nov 6 '13 at 12:06