clear;
nrSamp=10000;
lam=1;
y=rand(1,nrSamp); %uniform distribution
x=-1/lam*(log(1-y)); %neg ex distribution
t=cumsum(x); %time of events
dt=5; %count interval
nrDt=t(nrSamp)/dt-mod(t(nrSamp)/dt,1); %number of intervals that can be completely filled
count=zeros(1,nrDt); % counters for all events up to one interval
for i=1:nrDt
    a=find(t<i*dt); %find index of last event that fits in the interval
    if isempty(a)
        count(i)=0;
    else    
        count(i)=a(end);
    end    
end
poiss=diff(count); %difference is what occurs in just one interval only
fprintf(1, 'mean: %.2f std: %.2f\n', mean(poiss), std(poiss));
plot(poiss);
eVec=0:max(poiss);
h=histc(poiss,eVec);
figure; bar(h);