% a simple poisson death process (decay)

clear
N = 50; %number of samples
lam = 1;  %decay rate (dimensionless, hence = 1)
t = zeros(1,N);
% inital data % number of initial particles
K0 = 40;
K = K0;
J = [0:K0];

j=2;
while K>0
    
x = rand(1,N);
t(j,:) = t(j-1,:)-log(x)/(lam*K);
K = K-1;
j=j+1;
end


T = sort(t(j-1,:));

n = [1:N];

pn = (1-exp(-lam*T)).^K0;



    s = 1;
 tt = reshape([t(:,s)';t(:,s)'],2*(K0+1),1);
  nn = reshape([ K0:-1:0; K0:-1:0],2*(K0+1),1);

 figure(1)
 plot(tt(2:end),nn(1:end-1,1),tt,K0*exp(-tt),'--','linewidth',2)
 xlabel('\lambda t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 title('Sample path','fontsize',20)
figure(2)
 semilogy(tt,nn,'linewidth',2)
 xlabel('\lambda t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 title('semilog plot: Particle number vs. time','fontsize',20)
  
 figure(3)
  
 hist(t(K0+1,:),50)
 xlabel('\lambda t','fontsize',20)
  
 title('Extinction time histogram','fontsize',20)
 
 figure(4)
dT = T(2:N)-T(1:N-1);
dp = pn(2:N)-pn(1:N-1);

% exact pdf:
pdf = lam*K0*(1-exp(-lam*T)).^(K0-1).*exp(-lam*T);
%
plot(T(1:N-1),dp./dT,T,pdf,'--','linewidth',2);
xlabel('time','fontsize',20)
ylabel('Approximate PDF','fontsize',20)
 figure(5)
plot(T,n/N,'*',T,pn,'k')%
xlabel('time','fontsize',20)
ylabel('Approximate CDF','fontsize',20)

