clf

// Battage de carte 

// Fonction de base 

function y=melange(x)
N=length(x) 
n=grand(1,1,'uin',1,N)
y=[x(2:n),x(1),x(n+1:N)]
endfunction

y=melange(1:10)


// Une trajectoire avec N=10 cartes de longueur T=10
N=10;
T=10;
y=1:N;
for t=1:T
    y(t+1,:)=melange(y(t,:));
end
y


// Plus difficile: pour avoir une représentation graphique convaincante,
// on va suivre la place de chaque carte d'ou la necessite de calculer l'inverse
// d'une permutation

function z=perminv(y)
N=length(y)
for i=1:N
    z(1,y(i))=i;
end
endfunction

z=perminv(y(1,:))
for i=2:(T+1),
    z=[z;perminv(y(i,:))];
end

// On crée le graphique permettant de suivre la trajectoire de chaque carte 



clf
for i=1:N
plot2d(0:T,z(:,i)')    
end
title('Position des cartes en fonction du temps')

// Temps de remontée de la carte N en position 1 - echantillon de taille K
N=52;
K=1000;
S=zeros(1,K);
timer();
for k=1:1000
y=1:N;
s=0;
while y(1)<N
    y=melange(y);
    s=s+1;
end
S(k)=s+1;
end
timer()
clf
histplot(100,S);
title('histogramme de S - taille d échantillon K=1000')

// Methode 2 avec les lois géométriques
N=52;
K=1000;
T=zeros(1,K);
for k=1:1000
t=1;
for i=1:N
    t=t+grand(1,1,'geom',1/N);
end
T(k)=t;
end
timer()
clf
histplot(100,T);
title('histogramme de T - taille d échantillon K=1000')

// Normalité asymptotique de T
clf
histplot(100,(T-mean(T))/sqrt(variance(T)));
x=-3:0.1:3
plot2d(x,1/sqrt(2*%pi)*exp(-x.^2/2))
title('Normalité asymptotique de T')
