How to simulate Poisson Distribution Process?

Question

Safia 2022-9-5

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process

评论： 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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2022-9-5

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1042305

在 MATLAB Online 中打开

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 个评论
显示 12更早的评论隐藏 12更早的评论

Star Strider 2022-9-5

编辑：Star Strider 2022-9-5

在 MATLAB Online 中打开

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.

.

Star Strider 2022-9-5

在 MATLAB Online 中打开

@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 2022-9-29

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

请先登录，再进行评论。

Answer 2

William Rose 2022-9-5

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1042325

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Safia 2022-9-5

thank you, but i have a single segment.

请先登录，再进行评论。

Answer 3

Image Analyst 2022-9-5

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1042340

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 个评论
显示 2更早的评论隐藏 2更早的评论

William Rose 2022-9-6

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

Safia 2022-9-6

@Image Analyst thank you.

请先登录，再进行评论。

Answer 4

Bruno Luong 2022-9-5

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1042415

编辑：Bruno Luong 2022-9-5

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

William Rose 2022-9-5

在 MATLAB Online 中打开

@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

fprintf('Mean=%.1f, variance=%.1f\n',mean(v),var(v));

Mean=3.1, variance=3.0

histogram(v)

Try it.

William Rose 2022-9-6

@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 2022-9-6

编辑：Bruno Luong 2022-9-6

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

Bruno Luong 2022-9-6

在 MATLAB Online 中打开

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 2022-9-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1063455

在 MATLAB Online 中打开

@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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Safia 2022-9-29

@Bruno Luong thank you very much!

请先登录，再进行评论。

Answer 6

Bruno Luong 2022-9-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795340-how-to-simulate-poisson-distribution-process#answer_1043420

编辑：Bruno Luong 2022-9-6

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Farshid R 2022-10-1

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,

请先登录，再进行评论。

How to simulate Poisson Distribution Process?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

14 个评论
显示 12更早的评论隐藏 12更早的评论

更多回答（5 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

How to simulate Poisson Distribution Process?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

14 个评论 显示 12更早的评论隐藏 12更早的评论

更多回答（5 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

5 个评论 显示 3更早的评论隐藏 3更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

14 个评论
显示 12更早的评论隐藏 12更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论