% Birth2.m simulates the 'pure birth process'. Starting with
% nstart individuals, each having probability p of splitting in each
% time step, we continue until population size reaches N=nnow. 
%
% A for-loop is used to run multiple simulations, all with the same
% value of starting population size (nstart) and birth probability
% p, and then the histogram of times to reach N=nnow is plotted.
% Also plotted are the sample paths for the simulations

% vector to hold times to reach nnow
!c:
etimes=zeros(100,1);

for irep=1:100 

  % parameters
  nstart=1; p=0.1;

  %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<100;
  % update the stochastic population	
	m=rand(nnow,1);
	nbirth=sum(m<p);
	nnow=nnow+nbirth;
	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);
disp([t nbarnow])
plot(tvals,nt)
hold on;
% plot(tvals,nbar,'+');
xlabel('Time'); ylabel('Population size'); title('Pure Birth Process');


end;
figure
hist(etimes);
xlabel('Time'); ylabel('Frequency'); title('Time to reach N=100');
mean(etimes), std(etimes)