// Modele Wright Fisher
//Martingale positive, bornee par N : converge ps, L2, L1...
//La limite Minfini est a valeurs dans (0,N)
//On a (par exemple par convergence dominee, ou par convergence L1)
//E[Minfini]=E[M0]=a
//Donc (1-P(A disparait))N=a et finalement P(A disparait)=1-a/N

// 3. Illustration du processus 
N=20
Xg=zeros(1,100);
Xg(1)=10;
for k=1:99
    Xg(k+1)=grand(1,'bin',N,Xg(k)/N);
end
plot2d(Xg)

//Illustration de la proba de disparition 

//
Pabs(10,20)//Estimation de la proba d'extinction pour a=10 et N=20
Pabs(50,200)//Estimation de la proba d'extinction pour a=50 et N=200
clf;
Pabs2(20);//Estimation de la proba d'extinction pour N=20 en fonction de a.

clf;
plot2d([0 10 20 30 40 50],[0 Pabs(10,50) Pabs(20,50) Pabs(30,50) Pabs(40,50),1],-1) //Pour N=50. On ne le fait pas pour toutes les valeurs de a...

//4. Version avec mutations
//Simulation du processus
clf
pAB=0.3;pBA=0.7;
N=20
Xg=zeros(1,1000);
Xg(1)=10;
for k=1:999
    Xg(k+1)=grand(1,'bin',Xg(k),1-pAB) + grand(1,'bin',N-Xg(k),pBA);
end
plot2d(Xg)
//Il n'y a plus de convergence presque sure (en terme de chaine de Markov, il n'y a plus d'etats absorbants)

muest=LoilimWF(20,pAB,pBA);
clf
plot2d([0:N],muest)
clf
N=60;
muest2=muestim(60,pAB,pBA);
plot2d([0:N],muest2)

pAB=0.5;pBA=0.5;
muest3=muestim(60,pAB,pBA);
plot2d([0:N],muest3,2)

pAB=0.9;pBA=0.1;
muest4=muestim(60,pAB,pBA);
plot2d([0:N],muest4,3)
