% Processus de Galton-Watson pour une loi de fécondité
% à support sur {0,1,2}
% Utilise la fonction probadis

p0=input('entrer la probabilité associée à la valeur 0 : ');
p1=input('entrer la probabilité associée à la valeur 1 : ');
p2=1-p1-p0;
N=input('entrer le nombre d''itérations N : ');
NR=input('entrer le nombre de réalisations NR : ');
m=0*p0+1*p1+2*p2;
if (m>1)
    q=p0/p2;
else
    q=1;
end
sprintf('la moyenne de la loi de fécondite est m=%f\n', m)
sprintf('la probabilité d''extinction de la population est q=%f\n', q)
x=[0:2];
p=[p0,p1,p2];
X=zeros(NR,N+1);
X(1:NR,1)=1;
for i=1:NR
    for j=1:N
        Y=0;
        for k=1:X(i,j)
            Y=Y+probadis(3,x,p);
        end
        X(i,j+1)=Y;
        if (X(i,j+1)==0)
            break;  % on sort de la boucle car la population 
                    % s'est éteinte a la j^e itération
        end
    end
end
plot(0:N,q*ones(N+1,1),'r-');
hold on
plot(0:N,sum(X==0)/NR) 
title('Processus de Galton-Watson pour une loi discrète');
xlabel('Nombre de générations');
ylabel('Probabilité d''extinction');
legend('Empirique','Theorique',4);
hold off
