%% Convergence presque sure vers une constante

X=rand(1,500);
M=zeros(1,500);
M(1)=X(1);
for i = 2:500
   M(i)=max(X(i),M(i-1));
end;
plot(M);
% Mn est croissante et converge ps vers 1

%% Illustration de la loi des grands nombres

U=rand(1,5000);
X=ceil(6*U);
Z=cumsum(X)./[1:5000];
plot(Z);
hold on
plot([1,5000],[3.5,3.5]),'r';
hold off
%%
n=[1:5000];
plot(n.^(0.1).*(Z-3.5));
figure
plot(n.*(Z-3.5));
figure
plot(sqrt(n).*(Z-3.5));
%%
U=rand(1,5000);
Z=cumsum(U)./[1:5000];
plot(Z);
hold on
plot([1,5000],[0.5,0.5]),'r';
hold off
%%
% Simulation de Cauchy par la methode d'inversion de la fonction de
% repartition
n=1000; X=rand(1,n);Y=-tan(pi*(X-0.5));
M=cumsum(Y)./[1:n];
plot(M)
%Numeriquement il n'y a pas convergence ps : la loi de Cauchy n'est pas L1
%En fait M suit toujours une loi de Cauchy
Mn=zeros(1,500);
for i =1:500
    n=1000; X=rand(1,n);Y=-tan(pi*(X-0.5));
    Mn(i)=mean(Y);
end;
Ms=sort(Mn);
plot(Ms(20:480),[20:480]/500);
hold on;
x=[-4:0.1:4];
plot(x,atan(x)/pi+1/2,'r');
hold off;
%% Convergence ps pas vers une constante
U=rand(1,100);
E=(U<1/2);
V=cumsum(E./2.^[1:100]);
plot(V);
hold on
U=rand(1,100);
E=(U<1/2);
V=cumsum(E./2.^[1:100]);
plot(V,'r');
% la limite est differente : elle est aleatoire
hold off
%% Illustration de la loi limite : on illustre une convergence en loi 
W=zeros(1,2000);
for i=1:2000
U=rand(1,100);
E=(U<1/2);
W(i)=sum(E./2.^[1:100]);
end
Ws=sort(W);
plot(Ws,[1:2000]/2000);
% la limite est uniforme sur [0,1]
%% Avec une autre loi pour les Ei
clf
U=rand(1,100);
E=(U<1/2)+(-1)*(U>1/2);
V=cumsum(E./2.^[1:100]);
plot(V);
hold on
U=rand(1,100);
E=(U<1/2)+(-1)*(U>1/2);
V=cumsum(E./2.^[1:100]);
plot(V,'r');
clf
for i=1:2000
U=rand(1,100);
E=(U<1/2)+(-1)*(U>1/2);
W(i)=sum(E./2.^[1:100]);
end
Ws=sort(W);
plot(Ws,[1:2000]/2000);
% la limite est uniforme sur [-1,1]
%% Loi uniforme sur le Cantor
clf
W=zeros(1,2000);
for i=1:2000
U=rand(1,100);
E=2*(U<1/2);
W(i)=sum(E./3.^[1:100]);
end
Ws=sort(W);
plot(Ws,[1:2000]/2000);
%% Polya
%% 3
clf;
R=Polya(500,1,1);
X=R./(2+[0:500]');
Y=1-X;
Z=X./Y;
hold on
plot(X,'r')
plot(Y)
plot(Z,'g')
% On observe une convergence presque sure de X,Y et Z, vers quelque chose
% d'aleatoire : il suffit de refaire l'experience.
%% 4 
clf;
X=Polyaloi(500,500,1,1);
Xs=sort(X);
plot(Xs,[1:500]/500); %on observe une convergence vers la loi uniforme si r=b=1
figure;
X=Polyaloi(500,500,4,6);
Xs=sort(X);
plot(Xs,[1:500]/500) % si on change r et b ce n'est plus vrai
%% 
% Rn est uniforme discrete sur {1,..,n+1} : cela se montre par recurrence
% en utilisant la formule des probabilites totales
% On peut l'illustrer par exemple pour n=4, avec un histogramme ou la fonction de
% repartition empirique :
clf;
R=zeros(3000,1);
for k=1:3000
    R(k)=1;
for l = 1:4
    R(k) = Urne(R(k),l-1,1,1); 
end
end
histo(R,5,1,1);
figure;
plotrep(R);
hold on;
stairs(1/5*[1:5],'r');