//*********************************************************
// TEXTE : ÉVOLUTION ALÉATOIRE DU GÉNOTYPE D'UNE POPULATION
// Scopos vol 11, p. 21.
//*********************************************************

function [P]=Matrice(N,p)
// Étant donné un entier N et un vecteur colonne p à 2N+1 composantes,
// retourne la matrice markovienne P (2N+1)x(2N+1) où
// P(i,j) = C_{2N}^j (p_i)^j (1-p_i)^{2N-j}

// N.B. : pour une loi binomiale comme celle-ci, l'utilisation
// de la fonction 'binomial' est *beaucoup* plus rapide,
// (voir texte sur l'endémie : 0endemie.sci).

// Calcul des coefficients binomiaux (C_2N^k) par triangle de Pascal :
B=[1];
for i=1:2*N, B=[B,0]+[0,B];end
// Construction de P (sans boucle).
// Principe : pour construire la matrice carrée n x n A(i,j)=i^j, il suffit
// de construire 2 matrices carrées n x n, I et J puis de faire I .^ J :
// L=ones(1,n); I = (1:n)' * L ; J = L' * (1:n); A = I .^ J;
// Pour calculer P, on pourrait faire
// I=p * ones(B);                  // matrice carrée des p_i
// J=ones(B') * (0:2*N);           // matrice carrée des exposants j
// I1=ones(I)-I;                   // matrice carrée des 1 - p_i
// J1=2*N*ones(J)-J;               // matrice carrée des exposants 2N-j
// B1=ones(B') * B;                // matrice carrée des coeff. bin.
// P= B1 .* (I .^ J) .* ( I1 .^ J1);
// Pour économiser la mémoire on ne stocke pas les matrices intermédiaires
// et on écrit (c'est plus efficace mais complètement illisible !) :
P= (ones(B') * B);
P= P .* ((p * ones(B)) .^ (ones(B') * (0:2*N))) ;
P= P .* ((ones(p*ones(B))-p*ones(B)) .^ (2*N*ones(ones(B')*(0:2*N))-ones(B')*(0:2*N)));
endfunction

// Calcul de l'espérance du temps d'absorption.
// On considère la matrice restreinte aux états transients
// PT=P(2:2*N,2:2*N);
// Le vecteur des durées moyennes d'absoption est solution
// du système matriciel (I - PT) m = (1) (cf. exo 7).
// Il y a 2 classes récurrentes, le SE propre relatif à la
// valeur propre 1 de P est de dim. 2 (cf. prop 46 du cours),
// ce qui permet d'affirmer que ce système a une solution unique :
// considérer le système  L * (I - P) = (0), on en connaît déjà
// deux solutions lin. indép. (a,0,...,0) et (0,...,0,b).
// Si il existait une solution L telle que L_i non nul pour
// un i entre 1 et 2N, ça ferait un EV de dim >= 3.
// m est donc l'unique solution du système 
// (eye(PT)-PT) * m = ones(2*N-1,1)
// La résolution s'écrit en scilab :
// m = (eye(PT)-PT) \ ones(2*N-1,1);
// N.B.: calculer la matrice inverse inv(eye(PT)-PT) 
// serait maladroit (cf. Gauss / Gauss Jordan) !

function [E]=geno1()
// Retourne une matrice à 2N-1 lignes et 2 colonnes
// donnant le quotient
// espérance de la durée d'extinction d'un des 2 gènes / taille population
//   col 1 : en partant d'une population équilibrée (au départ
//           (N gènes A et N gènes a)
//   col 2 : en partant d'une population où les gènes A sont largement
//           minoritaires (x=5% de la population initiale)
// La fonction affiche également le graphe de ces 2 quotients
// en fonction de N (d'après l'article ces quotients devraient
// être à peu près constants...)
// N.B.: Pour N=194 on a un dépassement de capacité sur lps.
//
x=0.05;               // pourcentage (minorité = 5% du total)
K=17;                 // nombre d'essais
d=10;                 // pas 
for N=d:d:K*d,        // Taille totale de la population = 2N
  p=(0:2*N)' ./ (2*N);// Vecteur colonne des p_i= i/2N (0 <= i <= 2N)
  P=Matrice(N,p);     // Matrice markovienne
  PT=P(2:2*N,2:2*N);  // Matrice restreinte aux états transients
  m = (eye(PT)-PT) \ ones(2*N-1,1); // résolution système (cf. ci-dessus)
  E(N/d,1)=m(N)/(2*N);  // quotient (cas équilibré)
  E(N/d,2)=m(ceil(x*2*N))/(2*N);// quotient (cas A minoritaire)
end;
xbasc();
xselect();
X1=[d:d:10*K]';
X=[X1,X1];
plot2d(X,E,[2,5]);
endfunction;

function [E,S,L]=geno2(N,u,v)
// Modèle avec mutations.
// Entrées : N (taille population = 2N), limite (N =200 environ)
//           u = proba de mutation a->A, IL FAUT u>0 !
//           v = proba de mutation A->a, IL FAUT v>0 !
// Retourne le vecteur colonne L, proba invariante de la chaîne (réc. irréd.),
// son espérance E et son écart-type S.
// L(i) = proba stationnaire pour qu'il y ait i gènes `a' (et 2N-i gènes `A').
// L est l'unique solution de L (I-P+A) = L A = (1,1,...,1) où A=ones(P). 
// Vecteur colonne des p_i :
p=(0:2*N)' .* (1-u) + (2*N:-1:0)' .* v;
p= p ./ (2*N) ;
P=Matrice(N,p);                  // Matrice markovienne
L=ones(p')/(eye(P)-P+ones(P));   // matrice ligne (proba invariante)
L=L';                            // Le même en colonne (meilleur affichage)
E=(0:2*N) * L;                   // Espérance
S=sqrt(((0:2*N) .^2) * L - E^2); // Écart-type
endfunction;

// Le nombre moyen de gènes `a' en régime stationnaire est
// E= 2N * v / (u+v).
// On peut le vérifier pour différentes valeurs des paramètres, ex. :
// N=100;u=0.01;v=0.02; 
// geno2(N,u,v)*(u+v)/(2*N*v)

// Vérification du résultat donné par geno2 dans le cas u=v=1/2.
// Le calcul donne que L est une loi Bin(2N,1/2) : L(i) = C_{2N}^i * 2^{-2N}.
// Jusqu'à N=25 environ l'écart est acceptable.
// Au delà, l'erreur absolue est faible (évidemment) mais l'erreur
// relative s'envole... (Err = 30.35 pour N=30 !!!).
function [Erel,Eabs]=verif2(N)
// Solution analytique :
////B(1)=cdfbin("PQ",0,2*N,.5,.5);
//for j=2:2*N+1,
//  B(j)=cdfbin("PQ",j-1,2*N,.5,.5)-cdfbin("PQ",j-2,2*N,.5,.5);
//end;
// ou bien (triangle de Pascal) :
B=[1]; for i=1:2*N, B=[B,0]+[0,B];end
B= B' ./ 4^N;
[E,S,L]=geno2(N,.5,.5);
Eabs=max(abs(B-L));
Erel=max(abs(B-L) ./ B);
endfunction