function cvps(n)
  X=grand(1,n,"nor",0,1)
  Y=exp(X.^2/2)
  M=[]
  for i=1:n
    M=-sort(-[M,X(i)])
    med(i)=M(floor(i)/2+1)
end
xbasc();
plot2d([1:n],med,5)
xtitle('Convergence de la mediane empirique')
endfunction

function cvps2(n)
  X=grand(1,n,"nor",0,1)
  Y=exp(X.^2/2)
  M=[]
  for i=1:n
    M=-sort(-[M,Y(i)])
    med(i)=M(floor(i)/2+1)
end
xbasc();
plot2d([1:n],med,5)
//plot2d([1:n],cumsum(Y)./[1:n],5)
xtitle('Convergence de la mediane empirique')
endfunction

/////////////////////////////////////////////
function medgauss(n,p)
M=zeros(1,p)
for i=1:p
  X=-sort(-grand(1,2*n-1,"nor",2,1))
  M(i)=X(1,n);
  Moy(i)=mean(X)
end
Y=[-4:.01:4]
xbasc()
histplot(50,sqrt(2*n-1)*(M-2))
plot2d(Y,exp(-(Y.^2)/%pi)/%pi,5)
xtitle('Convergence de la mediane empirique pour des gaussiennes')
varmed=st_deviation(M)^2
varmoy=st_deviation(Moy)^2
disp('Variance empirique pour la mediane : '+string(varmed))
disp('Variance empirique pour la moyenne : '+string(varmoy))
endfunction


/////////////////////////////////////////////
function medlapl(n,p)
M=zeros(1,p)
Moy=zeros(1,p)
for i=1:p
  X=grand(1,2*n-1,"exp",0.5).*(2*(rand(1,2*n-1)>.5)-1)
  X=-sort(-X)
  M(i)=X(1,n);
  Moy(i)=mean(X)
end
Y=[-4:.01:4]
xbasc()
histplot(50,2*sqrt(2*n-1)*M)
plot2d(Y,exp(-Y.^2/2)/(sqrt(2*%pi)),5)
xtitle('Convergence de la mediane empirique pour des doubles exponentielles')
varmed=st_deviation(M)^2
varmoy=st_deviation(Moy)^2
disp('Variance empirique pour la mediane : '+string(varmed))
disp('Variance empirique pour la moyenne : '+string(varmoy))
endfunction

// Illustration du thm de Kolmogorov-Smirnov sur les fonctions de repartition
function kolmo(n,N)
D=zeros(1,N)
for i=1:N 
  X=-sort(-rand(1,n))
  D(i)=max(max(abs([1/n:1/n:1]-X)),max(abs([0:1/n:1-1/n]-X)))
end
D=-sort(-D)
xbasc()
plot2d2(sqrt(n)*[D(1) D],[0:1/N:1],5)
Y=[0.1:.01:sqrt(n)*D($)]
Yd=Y.^2
f=ones(Y)
for i=1:20000
  f=f+((-1)^i)*2*exp(-2*i^2*Yd)
end
plot2d(Y,f,2)  
xtitle('Thm de Kolmogorov-Smirnov')
legends(['F de rep KS','F de rep empirique'],[2,5],2)
disp(f((Y==1.36)))
endfunction

// Puissance du test de K-S pour des expo
function test(N,n)
eps=[0:.05:.7]
leps=length(eps)
c=zeros(eps)
for j=1:leps
  for i=1:N 
    X=-sort(-grand(1,n,"exp",1+eps(j)))
    FX=1-exp(-X)
    D=max(max(abs([1/n:1/n:1]-FX)),max(abs([0:1/n:1-1/n]-FX)))
    c(j)=c(j)+(sqrt(n)*D>1.36)
  end  
  c(j)=c(j)/N
  disp('Proportion de temps restant '+string(1-j/leps))
end
xbasc()
plot2d(eps,c,5)
xtitle('Puissance de Kolmo')
endfunction


// Comparaison des puissances de KS et CVM avec expo ou mélanges d'expo
function tests(N,n)
eps=[0:.05:.7]
leps=length(eps)
c=zeros(eps)
d=zeros(eps)
for j=1:leps
  for i=1:N 
    X=grand(1,n,"exp",1+eps(j))
//    X=grand(1,n,"exp",1)./(1+1*(rand(1,n)<eps(j)))
    X=-sort(-X)
    FX=1-exp(-X)
    D=max(max(abs([1/n:1/n:1]-FX)),max(abs([0:1/n:1-1/n]-FX)))
    c(j)=c(j)+(sqrt(n)*D>1.36)
    I=1/(12*n)+sum(((2*[1:n]-1)/(2*n)-FX).^2)
    d(j)=d(j)+(I>0.46)
  end  
  c(j)=c(j)/N
  d(j)=d(j)/N
  disp('Proportion de temps restant '+string(1-j/leps))
end
xbasc()
plot2d(eps,c,5)
plot2d(eps,d,2)
xtitle('Puissances comparées de Kolmo et Cramer')
legends(['Kolmo','Cramer'],[5,2],2)
endfunction



//////////////////////////////////////////
function pontbrownien(n)
U=rand(1,n); U=[U U]
U=-sort(-U);U=[0 U 1];
V=[0:1/n:1 ];V=[V V];V=-sort(-V)
DIF=sqrt(n)*(V-U);
xbasc();
plot2d(U,DIF,5)
plot2d([0,1],[0 0])
endfunction




//////////////////////////////////////////
function pont(U,i)
n=length(U)
U=rand(1,n); U=[U U]
U=-sort(-U);U=[0 U 1];
V=[0:1/n:1 ];V=[V V];V=-sort(-V)
DIF=sqrt(n)*(V-U);
plot2d(U,DIF,i)
plot2d([0,1],[0 0])
endfunction


////////////////////////////////////////
function verspont
xbasc()
U=rand(1,1000)
pont(U(1:10),3)
pont(U(1:100),6)
pont(U,5)
legends(['n=10','n=100','n=1000'],[3,6,5],1)
xtitle('Vers le pont brownien')
endfunction


/////////////////////////////////////
function test1(n,a)
m=2
X=grand(1,n,"exp",m)
tc=1/mean(X);
X=-sort(-X);
ma1=max(abs([(1/n):(1/n):1]-(1-exp(-tc*X))))
ma2=max(abs([0:(1/n):1-1/n]-(1-exp(-tc*X))))
D=sqrt(n)*max(ma1,ma2);
if D<=a then disp('Hypothese retenue');
  else disp('Hypothese rejetee') 
end
endfunction



/////////////////////////////////////
function test2(n,b)
m=2
X=grand(1,n,"exp",m)
S=cumsum(X)/sum(X);
S=S(1:$-1);
ma1=max(abs([(1/(n-1)):(1/(n-1)):1]-S));
ma2=max(abs([0:(1/(n-1)):1-1/(n-1)]-S));
D=sqrt(n)*max(ma1,ma2);
if D<=b then disp('Hypothese retenue');
  else disp('Hypothese rejetee') 
end
endfunction



/////////////////////////////////////
function niveau(a,b,p)
n=10;
m=2
c1=0;
c2=0;
for i=1:p
//tirage de lechantillon
  X=grand(1,n,"exp",m)
// test1
  tc=1/mean(X);
  T=-sort(-X);
  ma1=max(abs([(1/n):(1/n):1]-(1-exp(-tc*T))))
  ma2=max(abs([0:(1/n):1-1/n]-(1-exp(-tc*T))))
  D=sqrt(n)*max(ma1,ma2);
  D=sqrt(n)*max(abs([(1/n):(1/n):1]-(1-exp(-tc*T))));
  if D<=a then c1=c1+1
  end    
// test2
  S=cumsum(X)/sum(X);
  S=S(1:$-1);
  ma1=max(abs([(1/(n-1)):(1/(n-1)):1]-S));
  ma2=max(abs([0:(1/(n-1)):1-1/(n-1)]-S));
  D=sqrt(n)*max(ma1,ma2);
  if D<=b then c2=c2+1  
  end
end
disp('niveau du test 1 : '+string(c1/p))
disp('niveau du test 2 : '+string(c2/p))
endfunction



/////////////////////////////////////
function puissance(p)
n=10
a=.92
b=1.21
m=2
c1=0;
c2=0;
v=[(1/(n-1)):(1/(n-1)):1]; v=[v 1]
for i=1:p
//tirage de lechantillon
//  X=grand(1,n,"chi",1)
  X=abs(grand(1,n,"nor",0,1))
// test1
  tc=1/mean(X);
  T=-sort(-X);
  ma1=max(abs([(1/n):(1/n):1]-(1-exp(-tc*T))))
  ma2=max(abs([0:(1/n):1-1/n]-(1-exp(-tc*T))))
  D=sqrt(n)*max(ma1,ma2);
  if D<=a then c1=c1+1
  end    
// test2
  S=cumsum(X)/sum(X);
  S=S(1:$-1);
  ma1=max(abs([(1/(n-1)):(1/(n-1)):1]-S));
  ma2=max(abs([0:(1/(n-1)):1-1/(n-1)]-S));
  D=sqrt(n)*max(ma1,ma2);
  if D<=b then c2=c2+1
  end
end
disp('puissance du test 1 : '+string(1-c1/p))
disp('puissance du test 2 : '+string(1-c2/p))
endfunction