%%% Processus de Galton-Watson avec loi de fécondité portée par {0,...,n}.

clear
p=[1/4,1/4,1/2];  % Loi de reproduction sur {0,1,2,...}
ng=20;            % Nombre de générations 
nr=100;           % Nombre de réalisations

n=length(p)-1;    % Nombre maximal d'enfants possibles
pr=[0,cumsum(p)]; 
m=sum([0:n].*p);  
si=sqrt(sum(([0:n].^2).*p)-m^2)
sprintf('Moyenne m = %f\n',m)
sprintf('Écart type sigma = %f\n',si)

Z=zeros(nr,ng+1);
Z(1:nr,1)=1;

for i=1:nr % réalisations
  for j=1:ng % générations
	U=rand(1,Z(i,j));
	X=zeros(1,Z(i,j));
	for k=1:Z(i,j) % simulation des enfantements de la génération j
	  for l=1:n+1  % Z(i,j) réalisations i.i.d. de loi P
		if ( (U(k)>=pr(l)) & (U(k)<pr(l+1)) )
		  X(k)=l-1;
		  break
		end
	  end % for l
	end % for k
	Z(i,j+1)=sum(X);	
	if Z(i,j+1)==0  
	  break % On sort de la boucle for j si la population s'éteint
	end
  end % for j
end % for i
