// Approximation binomiale par Poisson

// Comparaison des deux lois
clf
Pb=binomial(0.01,100);
Ppo=exp(-1) ./factorial([0:20]);
plot2d3([0:20],Pb(1:21));
plot2d3([0.05:20.05],Ppo,3); // On decale legerement les abscisses pour voir quelque chose

// n trop petit
clf
Pb=binomial(0.1,10);
Ppo=exp(-1) ./factorial([0:10]);
plot2d3([0:10],Pb);
plot2d3([0.05:10.05],Ppo,2);


// p trop grand
clf
Pb=binomial(0.3,100);
Ppo=exp(-30)*30 .^[0:100] ./factorial([0:100]);
plot2d3([0:100],Pb);
plot2d3([0.05:100.05],Ppo,2);


// Comparaison des fonctions de repartition theorique
clf
Fb=cumsum(binomial(0.01,100));
Fpo=cumsum(exp(-1)./factorial([0:20]));
plot2d2([0:20],Fb(1:21));
plot2d2([0:20],Fpo,2);


// Illustration par fonction de repartition empirique
clf
Xb=zeros(1,10000);
for k = 1:10000
    Xb(k)=sum(rand(1,100)<0.01);
end;
Xbr=-gsort(-Xb);
plot2d2(Xbr,[1:10000]/10000);
Ppo=exp(-1)./factorial([0:10]);
plot2d2([0:10],cumsum(Ppo),2);



// Illustration par fonctions de repartitions empiriques si on ne connai pas la loi limite
clf
Xb=zeros(1,10000);
for k = 1:10000
    Xb(k)=sum(rand(1,100)<0.01);
end;
Xbr=-gsort(-Xb);
plot2d2(Xbr,[1:10000]/10000);
Xb2=zeros(1,10000);
for k = 1:10000
    Xb2(k)=sum(rand(1,500)<1/500);
end;
Xbr2=-gsort(-Xb2);
plot2d2(Xbr2,[1:10000]/10000,2);

// Approximation Binomiale par loi normale

// Comparaison des deux fonctions de repartition
clf
t=[-3:0.1:3];
n=100; 
p=0.4;
plot2d2(t,pbinom(t*sqrt(n*p*(1-p))+n*p,n,p));
plot2d(t,pnorm(t),2);


// n trop petit
clf

n=20; 
p=0.2;
a=-sqrt(n*p)/sqrt(1-p);
b=sqrt(n*(1-p))/sqrt(p);
t=[a:0.1:b];
plot2d2(t,pbinom(t*sqrt(n*p*(1-p))+n*p,n,p));
plot2d(t,pnorm(t),2);


// p trop petit

// Avec un echantillon
clf
Xb=zeros(1,10000);
n=1000; p=0.3;
for k = 1:10000
    Xb(k)=sum(rand(1,n)<p);
end;
Z=(Xb-n*p)/(sqrt(n*p*(1-p)));
Zr=-gsort(-Z);
plot2d2(Zr,[1:10000]/10000);
plot2d(t,pnorm(t),2);



// On peut aussi regarder la loi d'une binomiale : 
// elle ressemble a la densite d'une gaussienne de meme esperance et variance
clf
//Version sans stixbox
Pb=binomial(0.2,50);
plot2d3([0:50],Pb);
u=[0:0.1:50];
m=0.2*50;
s2=50*0.2*0.8;
plot2d(u,1/sqrt(2*%pi*s2)*exp(-(u-m)^2 ./(2*s2)),3);


//Version avec stixbox
clf;
Pb=dbinom([0:50],50,0.2);
plot2d3([0:50],Pb);
plot2d([0:0.1:50],dnorm([0:0.1:50],50*0.2,sqrt(50*0.2*0.8)),3)


// Si n est trop petit
//sans stixbox
clf;
Pb=binomial(0.2,20);
plot2d3([0:20],Pb);
u=[0:0.1:20];
m=0.2*20;
s2=20*0.2*0.8;
plot2d(u,1/sqrt(2*%pi*s2)*exp(-(u-m)^2 ./(2*s2)),3);

//avec stixbox
clf
Pb=dbinom([0:20],20,0.2);
plot2d3([0:20],Pb);
plot2d([0:0.1:20],dnorm([0:0.1:20],20*0.2,sqrt(20*0.2*0.8)),3)


// Si p est trop grand
//sans stixbox
clf;
Pb=binomial(0.99,50);
plot2d3([0:50],Pb);
u=[0:0.1:50];
m=0.99*50;
s2=50*0.99*0.01;
plot2d(u,1/sqrt(2*%pi*s2)*exp(-(u-m)^2 ./(2*s2)),3);

//avec stixbox

clf
Pb=dbinom([0:50],50,0.99);
plot2d3([0:50],Pb);
plot2d([0:0.1:50],dnorm([0:0.1:50],50*0.99,sqrt(50*0.99*0.01)),3)


// D'autres convergences

// Max d'exponentielles
clf
Y=zeros(1,1000);
for k=1:1000   
X=-log(rand(1,500))/2;
Y(k)=2*max(X)-log(500);
end;
plot2d2(-gsort(-Y),[1:1000]/1000);
t=[-5:0.1:5];
plot2d(t,exp(-exp(-t)),2);


// Min d'uniformes
clf
Y=zeros(1,1000);
for k=1:1000   
X=rand(1,500);
Y(k)=500*min(X);
end;
plot2d2(-gsort(-Y),[1:1000]/1000);
t=[0:0.1:5];
plot(t,1-exp(-t),3);


// Theoreme central limite
clf
Y=zeros(1,1000);
for k =1:1000
    X=rand(1,500);
   Y(k)=(mean(X)-0.5)*sqrt(500)*sqrt(12);
end;
plot2d2(-gsort(-Y),[1:1000]/1000);
plot2d([-3:0.1:3],pnorm([-3:0.1:3]),3);


// Avec un de
clf
Y=zeros(1,1000);
s=91/6-49/4;
for k =1:1000
    X=ceil(rand(1,500)*6);
   Y(k)=(mean(X)-3.5)*sqrt(500)/sqrt(s);
end;
plot2d2(-gsort(-Y),[1:1000]/1000);

plot2d([-3:0.1:3],pnorm([-3:0.1:3]),3);


// Avec Pareto
clf
Y1=zeros(1,1000);
for k=1:1000
X=rpareto(800,1.5);
Y1(k)=sqrt(800)*(mean(X)-1/3);
end
plot2d2(-gsort(-Y1),[1:1000]/1000);
Y2=zeros(1,1000);
for k=1:1000
X=rpareto(5000,1.5);
Y2(k)=sqrt(5000)*(mean(X)-1/3);
end
plot2d2(-gsort(-Y2),[1:1000]/1000,3);

//Avec une chaine de Markov

clf
d=5;
X1000=zeros(1,500);
X1100=zeros(1,500);
for k=1:500
    Y=2;
    for j=1:1000
        Y=Ehrenfest(d,Y);
    end;
    X1000(k)=Y;
    for j=1:100
        Y=Ehrenfest(d,Y);
    end;
    X1100(k)=Y;
end;
plot2d2(-gsort(-X1000),[1:500]/500);
plot2d2(-gsort(-X1100),[1:500]/500,3);
plot2d2([0:5],cumsum(binomial(0.5,5)),4)
