Poisson random number generator

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

Answer 1

[Andreas Goser]

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.

[Zhuofan Zheng]

great thanks

Answer 2

[Derek O'Connor]

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]

thank you for your help....

this one is quite helpful...

Answer 3

[Derek O'Connor]

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.  [e-mail]
%
  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

Answer 4

[Derek O'Connor]

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.  [e-mail]
%
  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]

an End is missing

[Ian Van Giesen]

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!

Answer 5

[Derek O'Connor]

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.  [e-mail]
%
  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]

Thank you for providing this, this is extremely helpful! Does anyone know why the poisrnd function is no longer in MatLab?

[John D'Errico]

Because the function name is spelled poissrnd, not poisrnd.

[Binlin Wu]

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.  [e-mail]
%
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

Answer 6

[Richard Willey]

Mark Steyvers has written a nice book titled "Computational Statistics with MATLAB" which can be downloaded from

http://psiexp.ss.uci.edu/research/teachingP205C/205C.pdf

The first chapter has a very good section describing inverse transform sampling which provides everything you need to know.

Answer 7

[Ahmed raheem]

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.

Answer 8

[Derek O'Connor]

@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.  [e-mail]
%

    k=1; produ = 1;
    produ = produ*rand;
    while produ >= exp(-lambda)
        produ = produ*rand;
        k = k+1;
    end
    X = k;

[Ahmed raheem]

thank you Derek.

it is performing well now...

i putted the code inside a (for loop) and the result was ok.

thanks again.

Answer 9

[Ahmed raheem]

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.  [e-mail]
%
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)

[Abdulramon Adeyiola]

@Ahmed raheem I tried using your code to generate 100 samples from Poisson distribution with parameter 1, but the mean of the resulting samples was way above 1. This is strange!

[Torsten]

L = exp(-lambda);
for i=1:n
    k=0; usave=1;
    while usave > L
        k = k+1;
        usave = usave*rand;
    end
    X(i) = k-1;
end