How to simulate Poisson Distribution Process?

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?

Answer 1

[Star Strider]

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

.

[William Rose]

@Safia,

@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)]);

    3.0167    2.7624

Good luck.

[Star Strider]

@William Rose — Thank you! (I very much wish that I was!)

[Safia]

Thank you @Star Strider and @William Rose for your reply, but the function 'poissrnd' requires Statistics and Machine Learning Toolbox. would you please suggest another solution

[Star Strider]

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)

v = 1×60

3 2 3 1 2 2 5 2 4 4 2 2 7 3 1 4 5 3 2 2 3 3 3 7 4 6 1 1 2 2

vm = buffer(v,10)

vm = 10×6

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

ixv = randperm(numel(vm))

ixv = 1×60

25 40 43 10 9 34 1 18 11 22 33 59 12 60 2 27 41 30 47 37 32 39 4 44 23 51 46 58 52 5

vnew = v(ixv)

vnew = 1×60

4 2 3 4 4 2 3 3 2 3 1 2 2 1 2 1 1 2 1 2 4 2 1 3 3 2 2 4 0 2

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 can't use this function "poissrnd", is there a free install of statistics and machine learning toolbox?

[Star Strider]

@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(:))

Values = 8×1

0 1 2 3 4 5 6 7

ixv = randperm(numel(vm));
vnew = vm(ixv)

vnew = 1×60

3 0 6 7 3 2 1 1 2 3 4 4 2 2 2 2 2 0 2 1 2 5 3 2 7 5 4 4 5 4

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]

many thanks @Star Strider , it works well.

[Star Strider]

As always, my pleasure!

[Safia]

@Star Strider Hi! i have installed new version of matlab therfore i found the statistics and machine learning toolbox and now i can use the function "poissrnd", i used the same code you provided but i want to keep the same sum of vnew.

before using the function,i just copied the vm ,so the sum still the same.

now i can't use copy of vm because the matrix is so large.

there is a new function that i can add?

thank you so much

lambda = 1/5;
v = poissrnd(lambda, 1, 3600);
vm = buffer(v,10);
ixv = randperm(numel(vm));
vnew = v(ixv);

[Star Strider]

@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]

@Star Strider thank you for reply! the problem is the vector v changed in each iteration , so the sum(v) usually changed. My question is there is another function to fix the sum because i want to use this value in matrix so i can't work with changeable values. Did you understand me?

[Star Strider]

I do not see any previous mention of sum(v) anywhere. What are you doing?

Would the poisscdf function be appropriate?

[Safia]

@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]

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.

Answer 2

[William Rose]

@Safia

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.

[Safia]

thank you, but i have a single segment.

Answer 3

[Image Analyst]

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.

[Safia]

it's just for one lane

[Image Analyst]

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]

@Image Analyst, thank you. Means a lot coming from you!

[Safia]

@Image Analyst thank you.

Answer 4

[Bruno Luong]

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)

[William Rose]

@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));

Mean=3.1, variance=3.0

histogram(v)

Try it.

[Bruno Luong]

@William Rose Just curious, what guide you to chose N=500?

[William Rose]

@Bruno Luong,

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]

@William Rose Thanks. No the answer is just perfect (N ~ lambda/0.01).

[Bruno Luong]

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)

Answer 5

[Bruno Luong]

@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

[Safia]

@Bruno Luong thank you very much!

Answer 6

[Bruno Luong]

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)

[Farshid R]

Hello @Bruno Luong,

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,