%branching process simulation

%offspring distribution on 0, 1, ... npmf

global npmf
global pmf

npmf=4;
pmf=zeros(1,npmf);
pmf(1)=0.3; pmf(2)=0.2; pmf(3)=0.1; pmf(4)=0.1; 

%simulation

ng=30; 
notrials=50;

for(tr=1:notrials)

pop=zeros(1,ng); 
for (i=1:ng)
    if (i==1) po=1; else po=pop(i-1);
    end
    u=rand(1,po);
    pmc=cumsum(pmf);
    pmcm=[0 pmc]; 
    z1=zeros(1,po);
    for (j=1:npmf)
        z1=z1+j*((u<=pmc(j))&(u>=pmcm(j)));
    end
    pop(i)=sum(z1);
    if ((notrials==1)&&(ng<=10))
        bar([1:po], z1);
        xt = [0];
        yt = [-0.5];
        str = sprintf('$X_{%d}=%d, X_{%d}=%d$', i-1, po, i, pop(i));
        text(xt,yt,str, 'Interpreter', 'Latex'); 
        ylim([-1,4]);
        title('Branching process: no. of offspring','FontSize',10);

        waitforbuttonpress;
        hold off;
    end
end

%plot([0:ng],[1 pop], 'LineWidth', 1);
mu=pmf*[1:npmf]';
plot([0:ng],[1 pop]./(mu.^[0:ng]), 'LineWidth', 1);
%ylim([0,100]);
hold on

end

xlabel('generation', 'Interpreter', 'Latex')
ylabel('population', 'Interpreter', 'Latex'); 
 
title('Branching process: population growth','FontSize',10);
 
hold off;
        
  