//Illustration de la convergence en loi via un histogramme

function convloi()
Z=[];
for i=1:1000
X=grand(1,1000,"def");  //loi uniforme sur [O,1[
Z=[Z (sqrt(1000*12)*(1/1000*sum(X)-1/2))];
end;
xbasc();
xtitle("Illustration du theoreme central-limit pour des variables aleatoires de loi U(0,1)");
histplot([-3:0.2:3],Z);
plot2d([-3:0.05:3],1/sqrt(2*%pi)*exp(-1/2*[-3:0.05:3].^2));
endfunction;

//Illustration de la loi des grands nombres

function lfgn(n,lambda)
X=grand(1,n,"exp",lambda);
xbasc();
xtitle("Illustration de la loi des grands nombres pour des variables aleatoires de loi E(lambda)");
plot2d([1:n],cumsum(X)./[1:n],3);
plot2d([1:n],lambda*ones(1,n),5);
endfunction;

//Exemple de simulation de variables aleatoires par methode du rejet

function rejet()
X=zeros(1,10000);
for i=1:10000;
x=rand(1,2);
while x(2) > sin(%pi*x(1)),
x=rand(1,2);
end,
X(i)=x(1);
end;
xbasc();
histplot(100,X);
Y=[0:0.05:1];
plot2d(Y,%pi*sin(%pi*Y)/2);
endfunction;


//Exemple de simulation de melange de lois

function melange(n)
X=grand(1,n,"exp",1/2);
Y=grand(1,n,"exp",1);
Z=grand(1,n,"def");
W=(Z<1/3).*Y+(Z>1/3).*X;
xbasc();
xtitle("Melange d une E(1/2) et d une E(1)");
histplot([0:0.05:5],W);
plot2d([0:0.05:5],1/3*exp(-[0:0.05:5])+4/3*exp(-2*[0:0.05:5]),5);
endfunction;




//Etude d'une chaine de Markov

//n=100 (longueur de la trajectoire) P=matrice de transition de la chaine X0=etat initial

function trajectoire(n,P,X0)
T=grand(n,"markov",P,X0);
X=[X0 T];
xbasc();
xtitle("Une trajectoire de la chaine de Markov");
plot2d(1:(n+1),X,5);
legends(["trajectoire"],[5],1);
endfunction;


//Pour une chaine irredusctible de noyau P, on illustre la convergence
// presque sure de la frequence empirique vers la valeur de la probabilite invariante
//X0 = etat initial  p1=valeur de pi(1)

function ergodique(n,P,X0,p1)
T=grand(n,"markov",P,X0);
X=[X0 T];
S=1.*(X==1);  //on etudie l'etat 1
V=cumsum(S);
xbasc();
xtitle("Illustration du theoreme ergodique");
plot2d(1:(n+1),V./(1:(n+1)),5);
plot2d(1:(n+1),p1*ones(1:(n+1)),2);
legends(["frequence empirique de 1";"p1"],[5;2],1);
endfunction;



//Etude de la convergence en loi de la chaine dans le cas ergodique
// b=mesure initiale

function convloi(n,P,b)
M=[b];
for j=2:n
M=[M ; M(j-1,:)*P];
end;
M=M';
xbasc();
xtitle("Convergence en loi de la chaine de Markov");
for k=1:8
plot2d(1:n,M(k,1:n),k);
end;
endfunction;

//Test du chi2 d'adequation de la loi de X_1 .. X_15 à la probabilite invariante (partant de X_0=1)
//Graphique donnant la p-valeur, et estimation de la probabilite de succes sous H0 par monte carlo
//chi2(15,MTF,1,u,100,100)

function chi2(N,P,X0,b,K,MC)
c=cdfchi("X",7,0.95,0.05);

U=[];
for n=1:N

res=[];

for k=1:MC



Z=[];

for i=1:K
S=grand(n,"markov",P,X0);
Z=[Z S(1,n)];
end;

R=[];
for j=1:8
u=sum(1*(Z==j));
R=[R u];
end;

D=sum(((K*b).^(-1)).*((R-K*b).^2));

res=[res 1*(D<c)];
end;

U=[U mean(res)];
end;

A=grand(1,200,"def");


xbasc();

subplot(211);
xtitle("Test du chi 2");
plot2d(1:N,U,5);
legends(["probabilite de succes sous H0"],[5],4);

subplot(212);
plot2d(1:200,((1:200).^(-1)).*cumsum(A),5);
plot2d(1:200,0.5*ones(1:200),1);
plot2d(1:200,((1:200).^(-1)).*cumsum(A)+1.96/sqrt(12)*sqrt((1:200).^(-1)),6);
plot2d(1:200,((1:200).^(-1)).*cumsum(A)-1.96/sqrt(12)*sqrt((1:200).^(-1)),6);

endfunction;


//Etude d'un processus de Poisson

// l := intensite du processus , t := intervalle de temps

function poisson(l,t)
N=grand(1,1,"poi",l*t);
T=grand(1,N,"unf",0,t);
M=-sort(-T);
xbasc();
xtitle("Trajectoire d un processus de Poisson");
plot2d2(M,(1:N),5);
legends(["nb de saut"],[5],4);
endfunction;

function res=poisson4(l,t)
N=grand(1,1,"poi",l*t);
T=grand(1,N,"unf",0,t);
M=-sort(-T);
res=M;
endfunction;


//etude asymptotique du processus de poisson

function asympt(l,T)
M=poisson4(l,T);
H=[1:T];
N=[];
for h=H
	N=[N,sum(M<=h)];
end;
xbasc();
xtitle("Etude asymptotique de Nt/t");
plot2d([0,H],[0,N./H],2);
plot2d2([0,H],l*ones([0,H]),5);
plot2d2(H,N./H-1.96*sqrt((N(T)/T)).*(H.^(-1/2)),6);
plot2d2(H,N./H+1.96*sqrt((N(T)/T)).*(H.^(-1/2)),6);
legends(["Nt/t";"lambda"],[2;5],4);
endfunction;


//etude de la loi de Nt a l'aide de la fonction de repartition empirique


function loint(l,t,MC)
D=[];

for i=1:MC
M=poisson4(l,t);
D=[D length(M)];
end;

S=-sort(-D);
X=[(S(1,1)-1) S (S($)+1)];
xbasc();
xtitle("Loi de Nt");
plot2d2([0:0.5:20],cdfpoi("PQ",[0:0.5:20],(l*t)*ones([0:0.5:20])),6);
plot2d2(X,[0:(1/(MC+1)):1],2);
legends(["f d r de la loi de poisson";"f d r empirique"],[6;2],4);
endfunction;


//Test du Chi2 d'independance

//realise un test d'independance , M la matrice des proportions
//On pourra reprendre les donnees de Tassi, Methodes statistiques

function res=chi2indep(M)
n=sum(M);
[a b]=size(M);
X=[];
Y=[];

for j=1:b
  Y=[Y sum(M(:,j))];
end;
for i=1:a
  X=[X; sum(M(i,:))];
end;

D=sum(((M - 1/n*(X*Y)).^2)./(1/n*(X*Y)));
d=(a-1)*(b-1);

k=cdfchi("X",d,0.95,0.05);

res=1*(D<k);

endfunction;



//Marche aleatoire et probleme de Dirichlet discret
//Voir Benaim, El Karoui, Promenade aleatoire

//sortie calcule le point de sortie d'une trajectoire dans M partant de (i0,j0)

function res=sortie(M,i0,j0)
i=i0;
j=j0;
while (M(i,j) < 2)
  X=grand(1,1,"def");
  if (X < 1/4) then 
    i=i+1; 
  elseif ((1/4 <= X)&(X < 1/2)) then 
    i=i-1;
  elseif ((1/2 <= X)&(X < 3/4)) then 
    j=j+1;
  else
    j=j-1;
   end,
end;
res=[i j];
endfunction;


//dirichlet resoud le probleme de Dirichlet dans M partant de (i0,j0) pour la fonction f, avec parametre MC

function res=dirichlet(M,f,i0,j0,MC)
Z=[];
for k=1:MC,
  X=sortie(M,i0,j0);  
  Z=[Z f(X(1),X(2))];
end;
res=mean(Z);
endfunction;

//resoud le probleme de Dirichlet discret dans le demi-plan superieur, compare le resulat avec une loi de Cauchy

function demiplan(N,MC)
R=[];
for i=1:MC
X=0;
Y=1;
while (Y>0)
  X=X+1/N*(2*(grand(1,1,"def") < 0.5)-1);
  Y=Y+1/N*(2*(grand(1,1,"def") < 0.5)-1);
end;
R=[R X];
end;
L=[-10:0.05:10];
xbasc();
xtitle("Etude de la loi de l abscisse de sortie");
plot2d(L,1/%pi*(1+L.^2).^(-1),5);
histplot([-10:0.5:10],R);
endfunction;

//resoud le probleme de Dirichlet discret sur le disque, compare le resultat avec celui donne par le noyau de Poisson

function sortie(x,y,r,MC,N) 
k=1; 
X=[x]; 
Y=[y]; 
while (((X($)^2+Y($)^2)<r^2) & (k<MC)) 
X=[X (X($)+1/N*(2*(grand(1,1,"def") <0.5)-1))]; 
Y=[Y (Y($)+1/N*(2*(grand(1,1,"def") <0.5)-1))]; 
k=k+1; 
end, 
L=linspace(- %pi,%pi,1000);
xbasc(); 
xtitle("Trajectoire d une marche aleatoire dans le disque");
plot2d(cos(r*L),sin(r*L),2); 
plot2d(X,Y,3); 
endfunction;

function a=disque(x,y,r,P,N,h,MC) 
A=[]; 
for i=1:MC 
k=1; 
X=[x]; 
Y=[y]; 
while (((X($)^2+Y($)^2)<r^2) & (k<N)) 
X=[X (X($)+1/P*(2*(grand(1,1,"def") <0.5)-1))]; 
Y=[Y (Y($)+1/P*(2*(grand(1,1,"def") <0.5)-1))]; 
k=k+1; 
end; 
A=[A h(X($),Y($))]; 
end; 
a=[1/MC*sum(A(1,:));1/MC*sum(A(2,:))]; 
endfunction;

function a=id(x,y) 
a=[x;y];
endfunction;

function y=noyaupoisson(r,theta,z) 
y=(r^2 ­ z(1)^2 ­ z(2)^2)./(z(1)^2 + z(2)^2 + r^2 + 2*r*sin(theta)*z(2)­ 2*r*z(1)*cos(theta)); 
endfunction;

function y=integralepoisson(z,r,N,f) 
for i=0:(N­1) 
l(i+1)=1/(N).*(f(r*cos(2*(%pi)*i/N), (r*sin(2*%pi*i/N))).*[noyaupoisson(r,2*%pi*i/N,z)]); 
end; 
y=sum(l);
endfunction;

function plotpoisson(r,z) 
theta=[0:0.05:2*%pi]; 
for i=1:max(size(theta)) 
N(i)=noyaupoisson(r,theta(i),z); 
end; 
plot2d(theta,N); 
endfunction;


//Convergence en loi, fonctions de répartition empiriques

//I un intervalle de valeurs pour le tracé

function plerf(I,m,v)
X=[];
for i=I
X=[X cdfnor("PQ",i,m,v)];
end;
plot2d(I,X,3);
endfunction;


//Compare la fonction de repartition empirique avec la fonction de repartition d'une gaussienne
//Trace un intervalle de confiance asymptotique fourni par le theoreme de Kolmogorov Smirnov

function frep(I)
xbasc();
xtitle("Convergence de la fonction de repartition empirique d'une gaussienne");
plerf((-3:0.05:3),0,1);
for i=I
	X=grand(1,i,"nor",0,1);
	Z=-sort(-X);
	plot2d2([-3 Z],(0:1/i:1));
	plot2d2([-3 Z],(0:1/i:1)+1.36/sqrt(i));
	plot2d2([-3 Z],(0:1/i:1)-1.36/sqrt(i));
end;
endfunction;

// Illustre le theoreme de Glivenko Cantelli

function glivenko(N)
X=grand(1,N,"nor",0,1);
G=[];
F=[];
for i=1:N
	F=[F cdfnor("PQ",X(1,i),0,1)];
end;	
for j=1:N
	S=-sort(-X(1,1:j));
	T=-sort(-F(1,1:j));
	G=[G max(max(abs(1/j*(1:j) - T),abs(1/j*(0:(j-1)) -T)))];
end;

xbasc();
xtitle("Illustration du theoreme de Glivenko Cantelli pour une gaussienne");
plot2d((1:N),G,5);
endfunction;


//Une etude de la convergence en loi pour Kolmogorov Smirnov

function kolmo(N,MC)
M=[];
for j=1:MC
X=grand(1,N,"exp",1)
S=-sort(-X);
F=1-exp(-S);
M=[M sqrt(N)*max(max(abs(1/N*(1:N)-F),abs(1/N*(0:(N-1))-F)))];
end;
xbasc();
xtitle("Illustration de la convergence en loi - Theoreme de Kolmogorov Smirnov");
histplot(20,M);
endfunction;