MATLAB Answers

## Poisson random number generator

Asked by Ahmed raheem

### Ahmed raheem (view profile)

on 6 Feb 2012
Latest activity Commented on by John D'Errico

### John D'Errico (view profile)

on 10 Jul 2019
Accepted Answer by Andreas Goser

### Andreas Goser (view profile)

Hi all please i need to know how to generate a Poisson distributed random variable without using the built-in function (poissrnd).

#### 0 Comments

Sign in to comment.

## 9 Answers ### Andreas Goser (view profile)

Answer by Andreas Goser

### Andreas Goser (view profile)

on 6 Feb 2012
Accepted Answer

If this is an acadamic exercise - you can look at the literature refererence
% References:
%  Devroye, L. (1986) Non-Uniform Random Variate Generation,
% Springer-Verlag.

Zhuofan Zheng

### Zhuofan Zheng (view profile)

on 26 Mar 2018
great thanks

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

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

Robert

### Robert (view profile)

on 14 Aug 2015
Thank you for providing this, this is extremely helpful! Does anyone know why the poisrnd function is no longer in MatLab?
John D'Errico

### John D'Errico (view profile)

on 10 Jul 2019
Because the function name is spelled poissrnd, not poisrnd.

Sign in to comment.

Answer by Richard Willey

### Richard Willey (view profile)

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.

#### 0 Comments

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

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.

Ahmed raheem

### Ahmed raheem (view profile)

on 7 Feb 2012
thank you for your help....
this one is quite helpful...

Sign in to comment.

Answer by Ahmed raheem

### Ahmed raheem (view profile)

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.

#### 0 Comments

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

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;

Ahmed raheem

### Ahmed raheem (view profile)

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.

Answer by Ahmed raheem

### Ahmed raheem (view profile)

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)

#### 0 Comments

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

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

#### 0 Comments

Sign in to comment.

Answer by Derek O'Connor

### Derek O'Connor (view profile)

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

PhD Student

### PhD Student (view profile)

on 10 Jul 2019
an End is missing

Sign in to comment.