% branch.m
% This MATLAB code simulates a single-type branching process
% and tracks the population size. You are prompted to supply
% the initial number in the population, the offspring 
% distribution as probabilities for producing
% 0, 1, 2, and 3 offspring (if this doesn't sum to 1 then 
% it is assumed that the reamining probability is assigned to 
% having 4 offspring), a time over which to
% run the process and the number of samples of the process 
% to make. 
%
% The output will be vectors: x contains a list of the
% final population size from each sample; xct1 contains
% a list of the final population size conditioned on
% non-extinction, that is for all samples which have not
% gone extinct. The extinction probability is also 
% produced, as well as a graph of all sample trajectories
%
!c:
n1=input('initial population size:    ')
p0=input('probability of 0 offspring:     ')
p1=input('probability of 1 offspring:     ')
p2=input('probability of 2 offspring:     ')
p3=input('probability of 3 offspring:     ')
n=input('Time period for process:     ')
q=input('how many simulations to run:     ')
p4=1-p0-p1-p2-p3;
%
% population size will be called pop
% and time will be ct. sumpop will contain 
% the pop size within each sample obtained by summing
% all the offspring. 
%
x=[ ];
etimes=[ ];
nt=[ ];
tvals=[ ];
for i=1:q
pop=n1;
et=n+1;
for ct=1:n
 sumpop=0;
  if pop==0
     sumpop=0;
     if et > ct 
          et=ct;
      else 
      end;    
  else for j=1:pop
      a=rand(1,1);
      if a >= p0+p1+p2+p3
        kids=4;
      elseif a >= p0+p1+p2
        kids=3;
      elseif a >= p0+p1
        kids=2;
      elseif a >= p0
        kids=1;
      else
        kids=0;
      end;
% end of the for j loop
     sumpop=sumpop+kids;
     end;
     pop=sumpop;
%
% store the results
%
        tvals=[tvals ct];
        nt=[nt pop];

   end;
end;
plot(tvals,nt),pause
hold on;
xlabel('Time'); ylabel('Population size'); title('Branching Process');
nt=[ ];
tvals=[ ];
etimes(i)=et;
x(i)=pop;
end;
hold off;
hist(x)
title('histogram for all samples at final time'),pause
hist(etimes)
title('histogram for extinction time'),pause
mean(etimes), std(etimes)