% Script file for Math 305 - lab assignment 7
% Simulation - Poisson process and Ch 7.2: 1, 5


% Simulate the 100 events satisfying a Poisson process with 
% rate lambda=2 events per unit of time.  Recall that the exponential cdf
% is p=F(t)= 1-e^(-\lambda t), 0\le t <1, and the inverse function is
% t=F^(-1)(t)= -(1/lambda))ln(1-p).  Let p be 100 uniform points from
% (0,1), then construct the corresponding t, and plot them.
%  Note here we use vectors instead of a "for loop".

lambda=2;
p=rand(1,100);
t=-(1/lambda).*log(1-p); % sequence of exponentially distr. times
ct=cumsum(t);            % sequence of partial sums of t's
y=zeros(1,100);          % want y coords to be 0
plot(ct,y,'*')       % each event time indicated by a *
% An alternative way of viewing these events, with the horizontal 
% axis time of even and vertical axis number of event:
stairs(ct,1:100)

% Now  count how many events in an interval of length T;
%  the average number should be lambda*T.  Here we use a "while loop"
% because we don't know in advance how many times the loop should run.
T=5;
tnow=0;
n=-1;
while tnow<T   % this will always be true the first time so n at least 0
p=rand();
t=-(1/lambda).*log(1-p);  
tnow=t+tnow;
n=n+1;    % the loop runs at least once so n >= 0
end

%  Now do 1000 such simulations 

lambda=2; T=5;
nlist=[];       %  A list for the number of t's in [0, T]
for i=1:1000
    tnow=0;
    n=-1;
    while tnow<T
        p=rand();
        t=-(1/lambda).*log(1-p);
        tnow=t+tnow;
        n=n+1;
    end
    nlist=[nlist,n];
end

mu=lambda*T;   % the mean number of events
b=2*mu+1;
x=0:b;
hist(nlist,x)  % histogram of numbers of events
% to compare to the Poisson graph we have to scale up
% the vertical axis of the pmf 

poi=exp(-mu)*(mu.^x)./factorial(x);  % Poisson pmf
hold on
 plot(x,1000*poi)
hold off
  
% See the handout (labwk7.pdf) for the rest of the assignment.