//****************************************
// PROCESSUS DE GALTON-WATSON À DEUX TYPES
// Texte Scopos vol 11, p. 95-102.
//****************************************

// Un modèle non-markovien créé à partir
// d'une chaîne de Galton-Watson à deux types.

// Paramètres :
q=0.5;   // proba de survie sans changement de type (type 1 ou 2)
p=1-q;   // proba de transformation (type 1 ou 2)
u=0.40;  // proba, pour un type 2, de transformation en 2 types 1
         // essayer u=0.24 (rho<1), u=0.25 (rho=1), u=0.26 (rho>1)
r=p-u;   // proba, pour un type 2, de mourir
if r<0,
  printf('*** Choisir u <= p ***');
  abort();
end
printf('\n p = %.2f, r = %.2f, u = %.2f , rho = %f', p, r, u, 1-p+sqrt(2*u*p));

// Nombre de générations étudiées
N=100;
// Taille de la génération initiale
U1=1;U2=0;    // 1 seul individu de type 1 au départ.

// Initialisations
t=([0:N])';   // vecteur colonne pour plot2d (temps)
X=zeros(t);   // vecteur colonne pour plot2d (taille génération t)
X(1)=U1+U2;

// Calcul du vecteur X = U1 + U2 :
for gen=1:N,
  // Pour chaque génération on tire 2 échantillons uniformes de tailles U1, U2
  // N.B. : le cas U1=0 ou U2=0 ne pose pas problème (R1=0 si U1=0).
  U=rand([1:U1],'uniform');
  R1=sum(bool2s(U<1-p));    // nombre de types 1 qui restent type 1
                            // U1-R1=nombre de types 1 transformés en type 2
  V=rand([1:U2],'uniform');
  R2=sum(bool2s(V<1-p));    // nombre de types 2 qui restent type 2
  T1=sum(bool2s(V>1-p+r));  // nombre de types 2 qui donnent 2 types 1
  U2=(U1-R1)+R2;
  U1=R1+2*T1;
  X(gen+1)=U1+U2;
end;

// Tracé de la trajectoire (X_n)
xbasc();         // effacer la fenêtre graphique
xselect();       // activer la fenêtre graphique (1er plan)
plot2d(t,X, style=2);