function bern
disp(1*(rand(1,10)<.5))

// v.a. de Poisson avec des expos
function poisson(n,l)
x=zeros(1,n)
for i=1:n 
  s=grand(1,1,'exp',1/l)
  c=0
  while s<=1 
    s=s+grand(1,1,'exp',1/l)
    c=c+1
  end
  x(i)=c
end  
emp=zeros(1,7)
for i=0:10 
  emp(i+1)=mean(x==i) 
end 
xbasc() 
plot2d2([0:10],emp,5)
plot2d2([0:10],exp(-l)*(l^[0:10])./([1 cumprod([1:10])]),2) 
xset("font size",3) 
xtitle('Simulation de variables aléatoires de Poisson grâce à des lois exponentielles')

// reflexion sur gaussiennes par rejet d'expo
function expgauss(c)
x=[-6:.01:6]
y=(c/2)*exp(-abs(x));
z=(exp(-(x.^2)/2))/sqrt(2*%pi)
xbasc();
plot2d(x,y,2)
plot2d(x,z,5)
//plot2d(x,y-z,3)

// loi des grands nombres unif
function lgn(n)
xbasc();
plot2d([0 n],[.5 .5],2)
plot2d([1:n],cumsum(rand(1,n))./[1:n],5)
xset("font size",3)
xtitle('LGN pour la loi uniforme sur [0,1]')

// loi des grands nombres pour 1/x^3
function lgnxt(n)
xbasc();
plot2d([0 n],[2 2],2)
plot2d([1:n],cumsum((sqrt(1-rand(1,n)))^(-1))./[1:n],5)
xset("font size",3)
xtitle('LGN pour une loi sans second moment')

// loi des grands nombres cauchy
function lgnc(n)
xbasc();
plot2d([0 n],[.5 .5],2)
plot2d([1:n],cumsum(tan(%pi*(rand(1,n)-.5)))./[1:n],5)
xset("font size",3)
xtitle('LGN pour la loi de Cauchy')

//tcl uniforme
function tcl(n,p)
z=zeros(1,p)
for i=1:p
  z(i)=sqrt(n)*(mean(rand(1,n))-.5)*sqrt(12)
end
xbasc();
histplot(30,z,2)
x=[-5:.01:5]
plot2d(x,(exp(-(x.^2)/(2)))/sqrt(2*%pi),5)
xtitle('TLC pour la loi uniforme sur [0,1]')

//tcl binomiale
function tclbin(n,p)
z=zeros(1,p)
for i=1:p
  z(i)=sqrt(n)*mean(2*(rand(1,n)<.5)-1)
end
xbasc();
histplot(30,z,2)
x=[-5:.01:5]
plot2d(x,(exp(-(x.^2)/(2)))/sqrt(2*%pi),5)
xset("font size",3)
xtitle('TLC pour la loi de Bernoulli')

//tcl expo
function tclexp(n,p)
z=zeros(1,p)
for i=1:p
  z(i)=sqrt(n)*(mean(grand(1,n,"exp",1))-1)
end
xbasc();
histplot(30,z,2)
x=[-5:.01:5]
plot2d(x,(exp(-(x.^2)/(2)))/sqrt(2*%pi),5)
xset("font size",3)
xtitle('TLC pour la loi exponentielle')

//tcl normale
function tclnor(n,p)
z=zeros(1,p)
for i=1:p
  z(i)=sqrt(n)*(mean(grand(1,n,"nor",0,1)))
end
xbasc();
histplot(30,z,2)
x=[-5:.01:5]
plot2d(x,(exp(-(x.^2)/(2)))/sqrt(2*%pi),5)
xset("font size",3)
xtitle('TLC pour la loi normale')

// nontcl pour 1/x^3
function nontcl(n,p)
z=zeros(1,p)
for i=1:p
  z(i)=sqrt(n)*(mean((sqrt(1-rand(1,n)))^(-1))-2)
end
xbasc();
histplot(30,z,2)
v=st_deviation(z)
x=[-10:.01:10]
plot2d(x,(exp(-(x.^2)/(2*v)))/sqrt(2*%pi*v),5)
disp('valeur maximale de l´échantillon : '+string(max(z)))
xset("font size",3)
xtitle('Mise en défaut du TLC  en l''absence de second moment')

// marche aleatoire sur Z
function marche(n)
xbasc();
plot2d([0:n],[0 cumsum(2*(rand(1,n)<.5)-1)],5)
plot2d([0 n],[0 0])
xset("font size",3)
xtitle('Trajectoire de la marche aléatoire simple')

// loi de l'arcsinus
function loiarcsin(n,p)
z=zeros(1,p)
for i=1:p
  x=cumsum(2*(rand(1,n)<.5)-1)
  z(i)=mean(1*(x>0)+.5*(x==0))
end
xbasc();
plot2d(-sort(-z),[1:p]/p,2)
x=[0:.01:1]
plot2d(x,asin(sqrt(x))/%pi,5)
xtitle('Loi de l´arcsinus')


// Tracé d'un échantillon de vecteurs gaussiens 
// de matrice de covariance donnée
function vecteur(n,K)
d=sqrt(length(K))
X=zeros(d,n)
A=chol(K)
for i=1:n
   X(:,i)=A'*grand(d,1,"nor",0,1)
end
xbasc()
plot2d(X(1,:),X(2,:),-1)


// melange de gaussienne
function melange(n,p,a,s)
y=2*a*(rand(1,n)<p)-a+grand(1,n,'nor',0,s)
//y=rand(1,n)
//y=[grand(1,sum(y<p),'nor',a,s),grand(1,sum(y>=p),'nor',-a,s)]
xbasc()
histplot(50,y,2)
x=[-6:.01:6]
plot2d(x,(p*exp(-((x-a).^2)/(2*s^2))+(1-p)*exp(-((x+a).^2)/(2*s^2)))/sqrt(2*%pi*s^2),5)
xset("font size",3)
xtitle('Distribution d''un mélange de deux gaussiennes avec p='+string(p)')


//estimation par moments
function esti(n,a,p)
y=2*a*(rand(1,n)<p)-a+grand(1,n,'nor',0,1)
ac=sqrt((cumsum(y.^2)./[1:n])-1)
pc=(ac+(cumsum(y)./[1:n]))./(2*ac)
xbasc();
plot2d([1:n],ac,2)
plot2d([1,n],[a,a],3)
plot2d([1:n],pc,5)
plot2d([1,n],[p,p],4)
legends(['Vrai a','Estimation de a','Vrai p','Estimation de p'],[3,2,4,5],4)
xset("font size",3)
xtitle('Estimation par la méthode des moments')