% TP5ex5 : Processus de Galton-Watson (avec loi binomiale)
% calcul de la probabilité d'extinction q

N=input('Nombre d''itérations : ');
n=2; % Premier paramètre de la loi binomiale.
p=input('Second paramètre de la loi binomiale : ');
NR=input('Nombre de réalisations : ');

m=n*p;
if m>1
   q=((1-p)/p)^2;
else
   q=1;
end
sprintf('Moyenne m=%f\n',m)
sprintf('Probabilité d''extinction q=%f\n',q)
X=zeros(NR,N+1);
X(1:NR,1)=1;
for i=1:NR
   for j=1:N
      X(i,j+1)=sum(rbinom([1,X(i,j)],n,p));
      if X(i,j+1)==0
         break % On sort de la seconde boucle for car la population s'éteint.
      end
   end
end
plot([0:N],sum(X==0)/NR)
hold on
plot([0:N],q*ones(N+1,1),'r-')
title('Processus de Galton-Walson Binomiale')
xlabel('Nombre de générations')
ylabel('Probabilités d''extinction')
grid
legend('Empirique','Théorique',4)
hold off

