How to simulate Poisson Distribution Process?
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?
Accepted Answer
14 Comments
@Star Strider is right, as always. The vector v has the number of cars passing in each minute. You can check that the sample mean and variance are approximately lambda, as expected.
lambda = 3;
v = poissrnd(lambda, 1, 60);
disp([mean(v),var(v)]);
Good luck.
Star Strider
on 5 Sep 2022
@William Rose — Thank you! (I very much wish that I was!)
Safia
on 5 Sep 2022
The probability density function of the Poisson Distribution is related to the exponential distribution. One approach would be to model that and then take its inverse.
An approach is described in the Wikipedia article on the Poisson distribution in the section on Random variate generation. You would have to adapt that code to MATLAB syntax.
Another alternative is to copy the result I posted and use it. If you want diffrerent versions (that is ‘shuffling’ it), use the randperm function on its indices:
lambda = 3;
v = poissrnd(lambda, 1, 60)
vm = buffer(v,10)
ixv = randperm(numel(vm))
vnew = v(ixv)
All the elements here are the same, only the order changes.
EDIT — (5 Sep 2022 at 17:17)
Discovered not all of ‘v’ was displaying, so used buffer to make all of it visible.
.
@Safia —
I was thinking that you could copy the ‘vm’ matrix and the use randperm to create randomised versions of it.
For example —
vm = [3 2 3 4 1 2
2 2 3 4 0 0
3 7 3 1 3 4
1 3 7 2 3 2
2 1 4 3 2 4
2 4 6 3 2 6
5 5 1 2 1 4
2 3 1 1 4 4
4 2 2 2 1 2
4 2 2 2 5 1];
Values = unique(vm(:))
ixv = randperm(numel(vm));
vnew = vm(ixv)
figure
histogram(vm(:), 8)
xt = xticks;
xticks(xt*8/9+1/2)
xticklabels(xt)
xlabel('Element Values')
ylabel('Number Of Values')
All I did here was to copy and paste ‘vm’ from my earlier code, then calculate ‘ixv’ here and create ‘vnew’ using it and the new randperm result. There are several repeats in ‘vm’ since the values go from 0 to 7, and there are 60 of them, so nothing is being missed.
.
Safia
on 6 Sep 2022
Star Strider
on 6 Sep 2022
As always, my pleasure!
Safia
on 29 Sep 2022
Star Strider
on 29 Sep 2022
@Safia — Now that you have the Statistics and Machine Learning Toolbox, you can create whatever values of Poisson random variables you want to with the poissrnd function, wherever and whenever you want to use them. My code was simply a work-around, since you needed the data and did not have other options.
Safia
on 29 Sep 2022
Star Strider
on 29 Sep 2022
I do not see any previous mention of sum(v) anywhere. What are you doing?
Safia
on 29 Sep 2022
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
2 votes
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
1 vote
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
Safia
on 5 Sep 2022
Image Analyst
on 5 Sep 2022
Yes, but you can do a million lanes to get a better, more precise distribution. Oh well I guess you don't understand. And I don't blame you because when I learned it in my "Probability, Statistical Optics, and Data Testing" course in grad school we spent a whole hour on this. But I think what @William Rose and I suggested is the way to go. His answer is probably better than mine. Good luck though.
William Rose
on 6 Sep 2022
@Image Analyst, thank you. Means a lot coming from you!
Safia
on 6 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
@Safia, going back to my eearlier suggesiotn about segmeents, I think you intrprtd "segment" diffeerently than I intended. A segmentin my suggestion is a short unit of time.
@Star Strider's approach is good - using the inverse of the known distribution function, taking uniform randoms numbersas looking up a random sample from the distribution.
My approach is more of a simulation. Examine each time short time segment and use a random number to decide whether or not one car passed by in that time window. This will yield a distribution of integers that approaches the Poisson distribtion, as the p=lambda/N becomes very small, where N=numbeer of timee segments per unit time. You then add up all the cars that pass (one at a time)in each unit of time.
lambda=3; %mean number of events per unit time
N=500; %number of segments in 1 unit of time
M=60; %number of time units to simulate
p=lambda/N; %event probability per segment of time
v=zeros(1,M);
for i=1:M
for j=1:N
if rand(1)<p,v(i)=v(i)+1; end
end
end
fprintf('Mean=%.1f, variance=%.1f\n',mean(v),var(v));
histogram(v)
Try it.
Bruno Luong
on 5 Sep 2022
Edited: Bruno Luong
on 5 Sep 2022
@William Rose Just curious, what guide you to chose N=500?
William Rose
on 6 Sep 2022
I chose N=500 because lambda=3, and I wanted p=lambda/N<0.01, and N=500 is the smallest "round number" that accomplishes that criterion. I wanted p<0.01 for the following reason: My simple simulation will only have 0 or 1 event per time step. It does not ever observe 2 or more events per time step. Therefore it is only an approimation. The chance of 2 events in a short segment is p^2/2; the chance of 3 events is p^3/3!, etc. (These numbers are very slightly higher than the actual probabilities of 2 and 3 events, but the diffeerence is only a factor of e^(-p), which is extremely close to 1.) By choosing p<0.01, we guarantee that the chance of two events in a segment (which my simulation fails to oberve) is < 5E-5, and the chance of 3 events in a segment is < 2E-7, etc. The low probability of missed multiple events means the simulaiton will be pretty accurate.
The nice thing about this way of simulating a Poisson process is that it reminds us of the simplicity of the underlying process. That simplicity is reflected in the fact that you only ever draw from the uniform distribution, with no calls to exponential functions or inverse CDFs, etc. The Poisson distribution arises naturally from this simple proceess, which is exact as one approaches the continuous limit of infinitesimally small time steps.
Sorry for overlong answer to the question.
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).
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)
@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,
Categories
Find more on Poisson Distribution in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
