How to simulate Poisson Distribution Process?

조회 수: 182(최근 30일)
Safia
Safia 2022년 9월 5일
댓글: Farshid R 2022년 10월 1일
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?

채택된 답변

Star Strider
Star Strider 2022년 9월 5일
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
Star Strider
Star Strider 2022년 9월 29일
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.

댓글을 달려면 로그인하십시오.

추가 답변(5개)

William Rose
William Rose 2022년 9월 5일
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.
  댓글 수: 1
Safia
Safia 2022년 9월 5일
thank you, but i have a single segment.

댓글을 달려면 로그인하십시오.


Image Analyst
Image Analyst 2022년 9월 5일
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
Safia
Safia 2022년 9월 6일
@Image Analyst thank you.

댓글을 달려면 로그인하십시오.


Bruno Luong
Bruno Luong 2022년 9월 5일
편집: Bruno Luong 2022년 9월 5일
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
Bruno Luong
Bruno Luong 2022년 9월 6일
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)

댓글을 달려면 로그인하십시오.


Bruno Luong
Bruno Luong 2022년 9월 29일
@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 2022년 9월 6일
편집: Bruno Luong 2022년 9월 6일
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
Farshid R
Farshid R 2022년 10월 1일
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,

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by