MATLAB Answers

Poisson random number generator

190 views (last 30 days)
Ahmed raheem
Ahmed raheem on 6 Feb 2012
Edited: Binlin Wu on 2 Apr 2021
Hi all please i need to know how to generate a Poisson distributed random variable without using the built-in function (poissrnd).

Accepted Answer

Andreas Goser
Andreas Goser on 6 Feb 2012
If this is an acadamic exercise - you can look at the literature refererence
% References:
% [1] Devroye, L. (1986) Non-Uniform Random Variate Generation,
% Springer-Verlag.

More Answers (8)

Derek O'Connor
Derek O'Connor on 6 Feb 2012
Dirk Kroese has excellent notes here: http://www.maths.uq.edu.au/~kroese/mccourse.pdf, which are based on his book:
D.P. Kroese, T. Taimre, Z.I. Botev: Handbook of Monte Carlo Methods. John Wiley & Sons, 2011.
His notes and book have lots of Matlab examples.
  1 Comment
Ahmed raheem
Ahmed raheem on 7 Feb 2012
thank you for your help....
this one is quite helpful...

Sign in to comment.


Derek O'Connor
Derek O'Connor on 7 Feb 2012
I prefer this:
% -------------------------------------------------------------
function S = PoissonSamp(lambda,ns);
% -------------------------------------------------------------
% Generate a random sample S of size ns from the (discrete)
% Poisson distribution with parameter lambda.
% Derek O'Connor, 6 Feb 2012. derekroconnor@eircom.net
%
S = zeros(ns,1);
for i = 1:ns
k=1; produ = 1;
produ = produ*rand;
while produ >= exp(-lambda)
produ = produ*rand;
k = k+1;
end
S(i) = k;
end

Derek O'Connor
Derek O'Connor on 3 May 2012
I would like to thank Kang Wook Lee of Berkeley for pointing out an error in the code above. The last line should be S(i) = k-1;
% -------------------------------------------------------------
function S = PoissonSamp(lambda,ns);
% -------------------------------------------------------------
% Generate a random sample S of size ns from the (discrete)
% Poisson distribution with parameter lambda.
% Fixed error: changed S(i) = k; to S(i) = k-1;
% Derek O'Connor, 3 May 2012. derekroconnor@eircom.net
%
S = zeros(ns,1);
for i = 1:ns
k=1; produ = 1;
produ = produ*rand;
while produ >= exp(-lambda)
produ = produ*rand;
k = k+1;
end
S(i) = k-1;
end
  2 Comments
Ian Van Giesen
Ian Van Giesen on 24 Jun 2020
Why was the last line changed? Was it to make the events where probability for 'success' a null event? Sorry in advance for a confusing question, still wrapping my head around this, but many thanks for the code!

Sign in to comment.


Derek O'Connor
Derek O'Connor on 24 Jul 2012
This is a cleaner fix of PoissonSamp
% -------------------------------------------------------------
function S = PoissonSamp3(lambda,ns);
% -------------------------------------------------------------
% Generate a random sample S of size ns from the (discrete)
% Poisson distribution with parameter lambda.
% Fixed error:
% CHANGED k = 1; produ = 1; produ = produ*rand
% TO k = 0; produ = rand;
% Derek O'Connor, 24 July 2012. derekroconnor@eircom.net
%
S = zeros(ns,1);
for i = 1:ns
k = 0;
produ = rand;
while produ >= exp(-lambda)
produ = produ*rand;
k = k+1;
end
S(i) = k;
end
  3 Comments
Binlin Wu
Binlin Wu on 1 Apr 2021
Would just like to make some minor changes to accept any array size:
function S = PoissonSamp3(lambda,varargin);
% -------------------------------------------------------------
% Generate a random sample S of size ns from the (discrete)
% Poisson distribution with parameter lambda.
% Fixed error:
% CHANGED k = 1; produ = 1; produ = produ*rand
% TO k = 0; produ = rand;
% Derek O'Connor, 24 July 2012. derekroconnor@eircom.net
%
nn = [varargin{:}];
S = zeros([nn,1]);
for i = 1:numel(S(:))
k = 0;
produ = rand;
while produ >= exp(-lambda)
produ = produ*rand;
k = k+1;
end
S(i) = k;
end

Sign in to comment.


Richard Willey
Richard Willey on 6 Feb 2012
Mark Steyvers has written a nice book titled "Computational Statistics with MATLAB" which can be downloaded from
The first chapter has a very good section describing inverse transform sampling which provides everything you need to know.

Ahmed raheem
Ahmed raheem on 7 Feb 2012
First of all, I'd like to thank you all for your cooperation with me. Based on Dirk Kroese, I wrote the following code to generate the Poisson random variable:
The m-file code:
n=1;
lambda=500;
for i=1:10000
x=rand(1);
a=1;
a=a*x;
if a>=exp(-lambda)
n=n+1;
continue
else
X(i)=n-1;
end
end
  • would you please check it if it's correct or not? i feel it is not correct.

Derek O'Connor
Derek O'Connor on 7 Feb 2012
@Ahmed, you're correct, it is not correct.
The function below is a Matlab translation of Kroese's algorithm. It seems to work ok but needs to be thoroughly tested. You should do this and let us know the results. Note that Poisson(L) ~ Norm(L,L), for large L.
The one thing I don't like about Kroese is his awful algorithm and programming style. Why does he use GOTOs instead of proper WHILEs etc?
% -------------------------------------------------------------
function X = Poisson(lambda);
% -------------------------------------------------------------
% Generate a random value from the (discrete) Poisson
% distribution with parameter lambda.
% Derek O'Connor, 6 Feb 2012. derekroconnor@eircom.net
%
k=1; produ = 1;
produ = produ*rand;
while produ >= exp(-lambda)
produ = produ*rand;
k = k+1;
end
X = k;
  1 Comment
Ahmed raheem
Ahmed raheem on 7 Feb 2012
thank you Derek.
it is performing well now...
i putted the code inside a (for loop) and the result was ok.
thanks again.

Sign in to comment.


Ahmed raheem
Ahmed raheem on 7 Feb 2012
This is the final code:
% -------------------------------------------------------------
function X = Poisson(lambda,n); % n represents the number of iterations
% -------------------------------------------------------------
% Generate a random value from the (discrete) Poisson
% distribution with parameter lambda.
% Derek O'Connor, 6 Feb 2012. derekroconnor@eircom.net
%
for i=1:n;
k=1; usave=1;
usave = usave*rand;
while usave >= exp(-lambda)
usave = usave*rand;
k = k+1;
end
X(i) = k;
end
hist(X)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!