// Shadock 
A=[5/6 1/12 1/12;1/4 1/2 1/4 ;1/4 0 3/4]
//Chaine finie irreductible aperiodique
//Donc une unique loi invariante, et convergence en loi

[vecp,valp]=spec(A');
muinv=vecp(:,1)'/sum(vecp(:,1))


mu0=[1,0,0]

mu0*A^5
mu0*A^10
mu0*A^50
mu0*A^100
//On observe bien une convergence en loi vers muinv

mu1=[1/2 1/4 1/4]
mu1*A^100 
//Encore convergence

A^5
A^10
A^100
//Convergence vers la matrice dont les lignes sont muinv


X=zeros(1,100);
X(1)=1;
for i=2:100
    X(i) = shadock(X(i-1));
end
plot2d2(X) //Sessin d'une trajectoire; 
//pour observer une convergence en loi, il en faut plusieurs.

Y=zeros(1,500);
for j=1:500
    Z=1;
        for i=2:100
        Z = shadock(Z);
        end
    Y(j)=Z;
end
loiYempirique=[sum(Y==1)/500,sum(Y==2)/500,sum(Y==3)/500]
clf
plot2d2([1:3],cumsum(loiYempirique));
plot2d2([1:3],cumsum(muinv),2)

// Zorro
//Chaine non irreductible, deux etats absorbants 0 et 8

A=zeros(9,9);
A(1,1)=1;
A(9,9)=1;
for i = 2:8
    A(i,i-1)=0.6;
    A(i,i+1)=0.4;
end
A

mu0=zeros(1,9);
mu0(4)=1;
mu0*A^3
mu0*A^10
mu0*A^100

mu=mu0*A^500;
mu(9) //Proba que Zorro soit sorti de prison au bout de 500 coups
// Donne une bonne Valeur approchee de proba que Zorro sorte de prison

//
N=0;
T=zeros(1,1000);
for k=1:1000
    X=3;
    T(k)=0;
while (X <8 & X>0)
    X=Zorro(X);
    T(k)=T(k)+1;
end
N=N+(X==8);
end
f=N/1000 //Estimation de la proba de sortir
s=sqrt(f*(1-f));
[f-s*1.96/sqrt(1000);f+s*1.96/sqrt(1000)] //Intervalle de confiance pour la proba de sortir
mean(T) //Estimation de l'esperance du temps total de jeu






// Ehrenfest
// 
A=MatEhrenfest(5)
//
Bin=dbinom([0:5],5,1/2)//utilise stixbox

Bin*A
// 3 La chaine est irreductible et de periode 2
// 4
X0=[1 0 0 0 0 0]
X20=X0*A^20
X50=X0*A^50
X51=X0*A^51
// La periode 2 fait que ces vecteurs ont des zeros soit aux places paires
// soit aux places impaires, donc il n'y a pas convergence vers la loi
// invariante. Mais on a bien convergence des soussuites X2n et  X(2n+1)

A^50
A^51
//La encore on voit que la suite A^n ne converge pas

X0=1/6*ones(1,6)
X20=X0*A^20
X50=X0*A^50
X51=X0*A^51
// Ici il y a bien convergence en loi vers la loi invariante

X=zeros(1,100);
for k=1:99
    X(k+1) = Ehren(X(k),5);
end
clf
plot2d(X)
//
Tv=0;
for k=1:1000
    Tv=Tv+vide(5);
end;
Tv=Tv/1000
1/Bin(1)
//
X=zeros(1,5000);
for k=1:4999
    X(k+1) = Ehren(X(k),5);
end
//m=0
clf
Erg0=cumsum(X==0)./[1:5000];
Bin0=Bin(1)*ones(1,5000);
plot2d(Erg0)
plot2d([1:5000],Bin0,2)
//m=1
clf
Erg1=cumsum(X==1)./[1:5000];
Bin1=Bin(2)*ones(1,5000);
plot2d(Erg1)
plot2d([1:5000],Bin1,2)
//m=2
clf
Erg2=cumsum(X==2)./[1:5000];
Bin2=Bin(3)*ones(1,5000);
plot2d(Erg2)
plot2d([1:5000],Bin2,2)
//m=3
clf
Erg3=cumsum(X==3)./[1:5000];
Bin3=Bin(4)*ones(1,5000);
plot2d(Erg3)
plot2d([1:5000],Bin3,2)
//m=4
clf
Erg4=cumsum(X==4)./[1:5000];
Bin4=Bin(5)*ones(1,5000);
plot2d(Erg4)
plot2d([1:5000],Bin4,2)
//m=5
Erg5=cumsum(X==5)./[1:5000];
Bin5=Bin(6)*ones(1,5000);
plot2d(Erg5)
plot2d([1:5000],Bin5,2)

//La chaine semble effectivement ergodique (voir theoreme...)

//B=Identite/(d+1) + d/(d+1)*A 
//B' et A' ont les memes vecteurs propres associes a la valeur propre 1
// donc les 2 modeles ont les memes lois invariantes
//La nouvelle chaine est aperiodique, donc cette foisci on a toujours
//convergence en loi vers la loi invariante


//Particules entrant et sortant d'un volume
clf
a=2;p=0.3;
X=Particules(p,a,1000,10);
plot2d(X)
k=[min(X):max(X)];
l=length(k);
prop=zeros(1,l);
for i=1:l
    prop(i)=sum(X==k(i))/1000;
end
clf
plot2d(k,prop,-1)
lam=a/p;
ppoi=zeros(1,l);
ppoi(1)=exp(-lam)*lam^(k(1))/factorial(k(1));
for i =2:l
    ppoi(i)=ppoi(i-1)*lam/(k(1)+i-1);//permet de calculer la loi du
    //Poisson avec peu de calculs.
end
plot2d3(k,ppoi,2)


// Modele genetique
//Deux etats absorbants : 0 et 2N
clf
N=10
Xg=zeros(1,100);
Xg(1)=N;
for k=1:99
    Xg(k+1)=grand(1,'bin',2*N,Xg(k)/(2*N));
end
plot2d(Xg)

clf
N=20
Xg=zeros(1,100);
Xg(1)=N;
for k=1:99
    Xg(k+1)=grand(1,'bin',2*N,Xg(k)/(2*N));
end
plot2d(Xg)

clf
N=100
Xg=zeros(1,1000);
Xg(1)=N;
for k=1:999
    Xg(k+1)=grand(1,'bin',2*N,Xg(k)/(2*N));
end
plot2d(Xg)


//
Pabs(10,20)
Pabs(20,40)
Pabs(100,200)
clf;
Pabs2(20);//Il semble raisonnable de poser Pabs=X0/(2N)

Pabs2(10); 
plot2d([0 10 20 30 40 50],[0 Pabs(10,50) Pabs(20,50) Pabs(30,50) Pabs(40,50),1],-1)

T=Tabs3();//Ca a l'air lineaire

clf;
Tabs2(3);
clf
Tabs2(10);
clf
Tabs2(200) //Ca a l'air d'etre une parabole

clf
u=0.3;v=0.7;
N=20
Xg=zeros(1,1000);
Xg(1)=N;
for k=1:999
    Xg(k+1)=grand(1,'bin',2*N,Xg(k)/(2*N)*(1-u)+(1-Xg(k)/(2*N))*(1-v));
end
plot2d(Xg)
//Plus d'etats absorbants
A=Matwright(10,u,v);
[vecp,valp]=spec(A');
muinv=vecp(:,1)'/sum(vecp(:,1))

A=Matwright(20,u,v);
[vecp,valp]=spec(A');
muinv=vecp(:,1)'/sum(vecp(:,1))

muest=muestim(20,u,v);
clf
plot2d([0:2*N],muinv)
plot2d([0:2*N],muest,2)
clf
N=60;
muest2=muestim(60,u,v);
plot2d([0:2*N],muest2)

u=0.5;v=0.5;
muest3=muestim(60,u,v);
plot2d([0:2*N],muest3,2)

u=0.9;v=0.1;
muest4=muestim(60,u,v);
plot2d([0:2*N],muest4,3)

//File d'attente
clf
X=zeros(1,1000);
X(1)=0;
for k=2:1000
    X(k)=Fattente(X(k-1),0.1,0.6);
end;
plot2d([1:1000],X);

X(1)=0;
for k=2:1000
    X(k)=Fattente(X(k-1),0.5,0.5);
end;
plot2d([1:1000],X,2);
X(1)=0;
for k=2:1000
    X(k)=Fattente(X(k-1),0.6,0.1);
end;
plot2d([1:1000],X,3);

// 
plotergo(0.1,0.6,10000);//la chaine semble raisonnablement ergodique
plotergo(0.5,0.55,10000);

//
convattente(0.1,0.6) //il y a bien convergence en loi vers la loi pi
convattente(0.5,0.5) // on n'observe pas de convergence en loi
convattente(0.6,0.1) // on n'observe pas de convergence en loi
