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

N=input('Nombre d''itérations : ');
p=input('Second paramètre de la loi géométrique : ');
NR=input('Nombre de réalisations : ');

m=(1-p)/p;
if m>1
   q=1/m;
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
   %i % pour connaitre lme degré d'avancement de la boucle
   for j=1:N
      X(i,j+1)=sum(rgeom([1,X(i,j)],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 géométrique')
xlabel('Nombre de générations')
ylabel('Probabilités d''extinction')
grid
legend('Empirique','Théorique',4)
hold off

