There is no difference between homogeneous and non-homogeneous Poisson distrib in this Matlab code

Question

I need to implement a traffic model with different probability distributions. I have this Matlab code (thanks to S. Dharmaraja) :

lambda=input('Enter The arrival Rate:');   % arrival rate
Tmax=input('Enter maximum time:');         % maximum time
clear T;
T(1)= 0;
a=input('Enter constant a>0:');
b=input('Enter constant b<1:');
S(1)=0;
i=1;

while T(i) < Tmax,
  U(i)=rand(1,1);
  T(i+1)=T(i)-(1/lambda)*(log(U(i))); % Homogeneous Poisson Distr.
  lambdat=a*(T(i+1))^(-b);
  u(i)=rand(1,1);
  if u(i)<=lambdat/lambda % <-- Here is my doubt/question
    i=i+1;
    S(i)=T(i);            % <-- Here is my doubt/question
  end
end
stairs(S(1:(i)), 0:(i-1));
title(['A Sample path of the non homogeneous Poisson process']);
xlabel(['Time interval']);
ylabel([' Number of event ']);

In this code, there are two lambdas. The first lambda "λ_rate" is the arrival rate, and the second lambda is λ(t) = a*t^-b.

As can be seen, there is a loop that is repeated until λ(t)/λ_rate >= rand(0,1). So what happens when this condition is met? S(i) = T(i), which means that the value obtained homogeneous poisson distribution is stored.

What is the difference between homogeneous and non-homogeneous? For me the difference is that in the second case there are more iterations than the first but in the end the graphs are equal.

What is your opinion?

PD; I put the code for the homogeneous case.

lambda=input('Enter The arrival Rate:');   % arrival rate
Tmax=input('Enter maximum time:');         % maximum time
clear T;
T(1)= 0;
i=1;

while T(i) < Tmax,
  U(i)=rand(1,1);
  T(i+1)=T(i)-(1/lambda)*(log(U(i)));
  i=i+1;
end

T(i)=Tmax;

stairs(T(1:(i)), 0:(i-1));
title(['A Sample path of the Poisson process with arrival rate ', num2str(lambda)])
xlabel(['Time'])
ylabel(['Number of Arrivals in [0,', num2str(Tmax), ']',])

Thanks,

Answer 1

I answered myself with the solution. This code was written by a student following the algorithm "2.6.3 Generating a nonhomogeneous Poisson Process" of the book "Simulation and Monte Carlo Method" of Rubinstein and Kroese (pg. 71).

The algorithm reads:

1. Set t=0, n=0 and i=0.
2. Increase i by 1.
3. Generate an independent random variable U(i) ~ U(0,1).
4. Set t = t - (1/λ)*ln(U(i)).
5. If t > T, stop; otherwise, continue.
6. Generate an independent random variable V(i) ~ V(0,1).
7. If V(i) <= λ(t)/λ, increase n by 1 and set T(n) = t. Go to step 2.

Where λ = max(λ(t)).

As λ(t) <= λ, I changed the function that brings the original code with the following:

lambdat = lambda*(1.2+cos(2*pi*T(i+1)/(Tmax/25)))/2.2;

where Tmax/25 is the period of the function.

The final code is:

lambda=input('Enter The arrival Rate:');   % arrival rate
Tmax=input('Enter maximum time:');         % maximum time
clear T;
T(1)= 0;
Per=input('Enter period of cosine: ');

T(1)=0;
n=0;
i=1;
while T(i) <= Tmax
    U(i) = rand(1,1);
    T(i+1) = T(i) - ((1/lambda)*log(U(i)));
    if T(i+1) > Tmax
        break;
    end
    lambdat = lambda*(1.2+cos(2*pi*T(i+1)/(Tmax/Per)))/2.2;

    V(i) = rand(1,1);
    if V(i) <= lambdat/lambda
        n = n+1;
        S(n) = T(i);
        Y(n)= lambdat;
    end
    i = i+1;
end
figure(1);
stairs(S(1:(n)), 0:(n-1));
title(['A Sample path of the non homogeneous Poisson process']);
xlabel(['Time interval']);
ylabel([' Number of event ']);

figure(2);
plot(S(:),Y(:),'*');
title(['lambda(t)']);

Answer 2

This is a technique known as "thinning," and was invented by Lewis & Shedler in 1978. You generate a homogeneous Poisson process via interarrival scheduling at rate lambda_max, the highest rate that occurs over the interval of interest. A given candidate event at time t is subjected to "thinning" by either accepting it with probability lambda(t) / lambda_max, or rejecting it and moving on to the next candidate. This adjusts the set of accepted event times to occur at the actual non-homogeneous arrival rates.

As an example, consider a non-homogeneous Poisson process with

            | 5 for t in [0,2),[3,5];
lambda(t) = |
            | 0 for t in [2,3).

With thinning you'll get an average of 10 occurrences for the first 2 time units, none at all between times 2 and 3, and then a resumption yielding an average of 10 for the remaining 2 time units. If you attempted to do this with instantaneous rates being used to generate the next arrival, you'd get one (forbidden) arrival somewhere between times 2 and 3, and then it wouldn't be able to restart later because with an instantaneous arrival rate of zero, your next arrival will occur in infinite time. (Okay, there's some chance of not getting an arrival between times 2 and 3, but the broader point stands that most of the time you'll get in trouble.)