How to simulate Poisson Distribution Process?

136 views (last 30 days)
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?

Accepted Answer

Star Strider
Star Strider on 5 Sep 2022
See if the poissrnd function will do what you want —
lambda = 3;
v = poissrnd(lambda, 1, 60)
v = 1×60
4 3 1 2 2 4 5 4 3 5 5 6 1 4 1 6 5 1 2 0 6 2 1 1 3 3 4 4 3 6
.
  14 Comments
Safia
Safia on 29 Sep 2022
@Star Strider sum(v) will be the sum of all cars rolling in the trafic during the period T.
v was a vector , each minute there are number of arrival cars. Now i want to fix the vector v, it means when i run many time the code, the v does not change or if values in v change, the sum still the same.
Star Strider
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.

Sign in to comment.

More Answers (5)

William Rose
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
Image Analyst on 5 Sep 2022
Do you want to simulate a stochastic process? Like let's say there are 1000 lanes of highway for the cars to travel on. These are the rows of a matric. And the columns are the times, from 1 to 60 (each column is a minute). And you populate the first column. Then in the next iteration (next minute) you move column 1 over to column 2 (because the cars in column 1 drove there) and now make a new column 1. Then in the next iteration you shift all the columns over one column and fill up column 1 with a few cars.
  4 Comments

Sign in to comment.


Bruno Luong
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
ans = 2.9968
std(r)^2 % ~ lambda
ans = 2.9934
histogram(r)
  5 Comments
Bruno Luong
Bruno Luong on 6 Sep 2022
Edited: Bruno Luong on 6 Sep 2022
@William Rose Thanks. No the answer is just perfect (N ~ lambda/0.01).
Bruno Luong
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
ans = 3.0034
std(r)^2 % ~ lambda
ans = 2.9986
histogram(r)

Sign in to comment.


Bruno Luong
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(:))
ans = 3
std(r(:))
ans = 3.0077

Bruno Luong
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)
ans = 3.0006
std(r)^2
ans = 3.0012
histogram(r)
  1 Comment
Farshid R
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,

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!