% DeathProcess2.m simulates the 'pure death process'. Starting with
% nstart individuals, each having probability p of dieing in each
% time step, we continue until all have died. 
%
% A for-loop is used to run multiple simulations, all with the same
% value of starting population size (nstart) and death probability
% p, and then the histogram of extinction times is plotted.

% vector to hold extinction times
etimes=zeros(100,1);

for irep=1:100 

  % parameters
  nstart=50; p=0.5;

  %initialize variables for current population size
  nnow=nstart; nbarnow=nstart; t=0;

  % initialize vectors storing the results over time
  nt=nstart; nbar=nstart; tvals=0;

  while nnow>0;
  % update the stochastic population	
	m=rand(nnow,1);
	ndie=sum(m<p);
	nnow=nnow-ndie;
	nbarnow=(1-p)*nbarnow;
   % store the results
   	t=t+1;
	tvals=[tvals;t];
	nt=[nt;nnow];
	nbar=[nbar;nbarnow];
  end; 

  etimes(irep)=t;
  disp(irep);

end;

hist(etimes);
mean(etimes), std(etimes)