//*********************************************
// PROCESSUS DE GALTON-WATSON AVEC IMMIGRATION
// Texte Scopos 1 p. 85.
//*********************************************

// Loi des Yn,i : elle ne charge que 0, 1 et 2.
p0=0.2;
p2=0.45;
p1=1-p0-p2;
m=p1+2*p2;
printf('\n Moyenne des Yn,i: m = %f\n', m);

// Loi des In : Bernoulli de paramètre 1/6
p=1/6;

// Nombre d'itérations:
N=40;

// Initialisation du processus
S=1;             // nombre initial de particules
X=zeros(N,1);    // vecteur des valeurs successives des X_n
X(1)=S;
rand('uniform');

// Construction du processus (X_n)
for i=2:N
  // Construction des Y_n,i:
  U=rand(1,S);
  Y=bool2s(U>p0) + bool2s(U>p0+p1);
  // Construction de I_n:
  V=rand(1);
  I=bool2s(V<p);
  // Construction de X_n:
  S=sum(Y)+I;
  // On stocke les valeurs des X_n dans le vecteur X
  X(i)=S;
end

// Calcul des (L_n)
L=X./(m.^[1:N]');

// Calcul de l'estimateur m_n de m : on n'a que n<N car v_n dépend de X_{n+1}
sN=cumsum(X(1:N-1));               // vecteur des s_n
sN1=cumsum(X(2:N));                // vecteur des s_{n+1}-1
tN=cumsum(X(1:N-1).^2);            // vecteur des t_n
vN=cumsum(X(1:(N-1)).*X(2:N));     // vecteur des v_n
detN=(N*(tN+1)-sN^2);              // vecteur des déterminants d_n
estN=(N.*vN - sN.*sN1) ./ detN;    // vecteur des estimateurs m_n

// Visualisation du comportement asymptotique de (X_n), (L_n) et (m_n) :
xbasc();
xselect();
// Fenêtre supérieure (Xn) en vert (3)
subplot(3,1,1);
plot2d([1:N]',X,style=3);
legends("Variations de X_n",style=3,opt="ul");
// Fenêtre inférieure (Ln) en bleu (2)
subplot(3,1,2);
plot2d([1:N]',L,style=2);
legends("Variations de L_n",style=2,opt="ul");

// Illustration de la convergence p.s. de m_n vers m.
subplot(3,1,3);
plot2d([1:N-1]', estN, style=5);
xsegs([1;N-1],[m;m], style=2);
legends("Convergence de m_n",style=5,opt="ul");