close all
T=10;      %minutes
lambda= 1; %customers per minute; 
nr_experiments=10^4;
n=1000;    

h=T/n;
p = lambda*h;

% Generate arrivals for all times and experiments
arrival = binornd(1,p,n,nr_experiments);

% Compare with Poisson pmfs
x=[0:30];
pdf_approx = hist(sum(arrival),x)/nr_experiments;
bar(x,pdf_approx)
hold on
plot(x,poisspdf(x,lambda*T),'r','Linewidth', 2)

figure
pdf_approx = hist(sum(arrival(1:n/2,:)),x)/nr_experiments;
bar(x,pdf_approx)
hold on
plot(x,poisspdf(x,lambda*T/2),'r','Linewidth', 2)

% Compute time of first arrival
time=0;
experiment=1;
time_histogram = zeros(n,1);
while experiment <= nr_experiments
    time = time+1;
    if arrival(time, experiment)
        time_histogram(time)=time_histogram(time)+1; 
        experiment = experiment+1;
        time=0;
    end
end

% Compare with exponential pdf
figure
plot([1:n]*h,time_histogram/nr_experiments/h)
hold on
plot([1:n]*h,exppdf([1:n]*h,lambda),'r')


figure
plot([1:n]*h,cumsum(time_histogram/nr_experiments))
hold on
plot([1:n]*h,expcdf([1:n]*h,lambda),'r')

