# How to simulate Poisson Distribution Process?

166 views (last 30 days)

Show older comments

Hello everyone!

I have a number of vehicles that will pass through a segment of highway, the arrival rate follows a Poisson distribution with lambda=3 (vehicle/min)

T=60 min, dt=1 min.

How can I simulate this process?

##### 0 Comments

### Accepted Answer

Star Strider
on 5 Sep 2022

##### 14 Comments

Star Strider
on 29 Sep 2022

I am not certain that I understand.

The Poisson process determines when the cars will arrive, not how many there will be. The requirement of 3 vehicles / minute means that the number of cars will be the same (or very nearly the same) for any specific measuring interval, providing the intervals are of the same length, and the length is, for example, several minutes.

### More Answers (5)

William Rose
on 5 Sep 2022

Yes. You can usee the binomial (or is it Bernoulli?) approximation: Chop each minute into N equal segments, where lambda/N <<1. In each segmeent, the chance of 1 car passing is p=lambda/N. The chance of two cars passing is negligible if N is large enough. That is why this is an approximation. So you "flip a coin" by drawing from rand() in each segment. Add 1 to the car count, each time rand() is less than p.

Then reinitialize your car counter to zero, and do it all again for the next minute.

The above is just an outline; I do not have time to check the details or providee code, but this will give a good approximation I think.

Try it and good luck.

Image Analyst
on 5 Sep 2022

Bruno Luong
on 5 Sep 2022

Edited: Bruno Luong
on 5 Sep 2022

If you don't have staistics toolbox

% Generate n random integers follow Poisson with parameter (==mean) lambda

lambda=3;

n = 1e6; % 60 in your case

% Find the max of possible value

k = 0;

tol = eps(exp(lambda));

p = 1;

while p > tol

k = k + 1;

p = p * lambda/k;

end

kmax = k;

% cumulative sum

K=0:kmax;

c=[0 cumsum(lambda.^K./factorial(K))]*exp(-lambda);

c=c/c(end);

r=discretize(rand(1,n),c)-1;

mean(r) % ~ lambda

std(r)^2 % ~ lambda

histogram(r)

##### 5 Comments

Bruno Luong
on 6 Sep 2022

The Taylor partial sum of the exponential can be computed by the incomplete gamma function:

% Generate n random integers follow Poisson with parameter (==mean) lambda

lambda=3;

n = 1e6; % 60 in your case

% Find the max of possible value

k = 0;

tol = eps(exp(lambda));

p = 1;

while p > tol

k = k + 1;

p = p * lambda/k;

end

kmax = k;

% cumulative sum

c = gammainc(lambda,0:kmax+1,'upper');

c = c/c(end); % make sure the last bin is 1.

r=discretize(rand(1,n),c)-1;

mean(r) % ~ lambda

std(r)^2 % ~ lambda

histogram(r)

Bruno Luong
on 29 Sep 2022

@Safia is you want integer with fixed sum, I propose this, and it is no longer Poisson law strictly speaking

lambda=3;

T=60; % period length, in minutes

n=1e3; % number of sequences of T

% r is n x T, random sum(r,2) is s

% so the "frequency" of r is s/T ~ lambda

s=round(lambda*T);

r=randi([0 s],n,T-1);

r=diff([zeros(n,1),sort(r,2),s+zeros(n,1)],1,2);

% r does not follow poisson, but not faraway as show here

histogram(r(:))

mean(r(:))

std(r(:))

Bruno Luong
on 6 Sep 2022

Edited: Bruno Luong
on 6 Sep 2022

Another way based on exponetial distribution, which is the time between 2 events

% Generate n random integers follow Poisson with parameter (==mean) lambda

lambda=3;

n = 1e6;

r = [];

tlast = 0;

while tlast < n+1

% we usualy need only one iteration

m = ceil((n+1-tlast)*lambda*1.1); % 10% of margin

s = cumsum(-log(1-rand(m,1))/lambda);

r = [r; tlast+s]; %#ok

tlast = r(end);

end

r = accumarray(ceil(r),1);

r(n+1:end) = []; % remove the exceed tail

% Check

mean(r)

std(r)^2

histogram(r)

##### 1 Comment

Farshid R
on 1 Oct 2022

I have an optimization problem(Consensus Optimization problem), unfortunately, I posted the problem on the MATLAB site,

@Sam Chak told me that I can send a message to you and I can discuss my problem with you.

I link to the question:

https://www.mathworks.com/matlabcentral/answers/1812615-optimization-with-fmincon-command-in-simulink?s_tid=mlc_ans_men_view&mentions=true#comment_2390935

In the question, I have given both simulink and relations completely and I discussed with @Sam Chak but it did not reach the desired result and @Sam Chak said to send you a message.

Please guide me . I am looking forward to your positive answer.

Best regards,

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!