%% Approximation binomiale par Poisson

%% Comparaison des deux lois
Pb=dbinom([0:20],100,0.01);
Ppo=exp(-1)./factorial([0:20]);
hold on
stem([0:20],Pb);
stem([0:20],Ppo,'r');
hold off
%% n trop petit
clf
Pb=dbinom([0:10],10,0.1);
Ppo=exp(-1)./factorial([0:10]);
hold on
stem([0:10],Pb);
stem([0:10],Ppo,'r');
hold off

%% p trop grand
clf
Pb=dbinom([0:100],100,0.3);
Ppo=exp(-30).*30.^[0:100]./factorial([0:100]);
hold on
stem([0:100],Pb);
stem([0:100],Ppo,'r');
hold off

%% Comparaison des fonctions de repartition theorique
clf
Fb=pbinom([0:20],100,0.01);
Fpo=cumsum(exp(-1)./factorial([0:20]));
hold on
stairs([0:20],Fb);
stairs([0:20],Fpo,'r');
hold off;

%% 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=sort(Xb);
stairs(Xbr,[1:10000]/10000);
hold on;
Ppo=exp(-1)./factorial([0:10]);
stairs([0:10],cumsum(Ppo),'r');
hold off;


%% Illustration par fonctions de repartitions empiriques si on ne connait
%% pas la loi limite
clf
Xb=zeros(1,10000);
for k = 1:10000
    Xb(k)=sum(rand(1,100)<0.01);
end;
Xbr=sort(Xb);
stairs(Xbr,[1:10000]/10000);
hold on;
Xb2=zeros(1,10000);
for k = 1:10000
    Xb2(k)=sum(rand(1,500)<1/500);
end;
Xbr2=sort(Xb2);
stairs(Xbr2,[1:10000]/10000,'r');

%% Approximation Binomiale par loi normale

%% Comparaison des deux fonctions de repartition
clf
t=[-3:0.1:3];
n=50; 
p=0.2;
stairs(t,pbinom(t*sqrt(n*p*(1-p))+n*p,n,p));
hold on;
plot(t,pnorm(t),'r');
hold off;

%% n trop petit
clf
t=[-3:0.1:3];
n=20; 
p=0.2;
stairs(t,pbinom(t*sqrt(n*p*(1-p))+n*p,n,p));
hold on;
plot(t,pnorm(t),'r');
hold off;

%% p trop petit
clf
t=[-3:0.1:3];
n=50; 
p=0.01;
stairs(t,pbinom(t*sqrt(n*p*(1-p))+n*p,n,p));
hold on;
plot(t,pnorm(t),'r');
hold off;

%% On peut aussi regarder la loi d'une binomiale : 
%% elle ressemble a la densite d'une gaussienne de meme esperance et variance
clf
Pb=dbinom([0:50],50,0.2);
hold on
stem([0:50],Pb);
plot([0:0.1:50],dnorm([0:0.1:50],50*0.2,sqrt(50*0.2*0.8)),'r')
hold off

%% Si n est trop petit
clf
Pb=dbinom([0:20],20,0.2);
hold on
stem([0:20],Pb);
plot([0:0.1:20],dnorm([0:0.1:20],20*0.2,sqrt(20*0.2*0.8)),'r')
hold off

%% Si p est trop grand
clf
Pb=dbinom([0:50],50,0.99);
hold on
stem([0:50],Pb);
plot([0:0.1:50],dnorm([0:0.1:50],50*0.99,sqrt(50*0.99*0.01)),'r')
hold off


%% 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;
stairs(sort(Y),[1:1000]/1000);
hold on
t=[-5:0.1:5];
plot(t,exp(-exp(-t)),'r');
hold off

%% Min d'uniformes
clf
Y=zeros(1,1000);
for k=1:1000   
X=rand(1,500);
Y(k)=500*min(X);
end;
stairs(sort(Y),[1:1000]/1000);
hold on
t=[0:0.1:5];
plot(t,1-exp(-t),'r');
hold off

%% 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;
stairs(sort(Y),[1:1000]/1000);
hold on
plot([-3:0.1:3],pnorm([-3:0.1:3]),'r');
hold off

%% 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;
stairs(sort(Y),[1:1000]/1000);
hold on
plot([-3:0.1:3],pnorm([-3:0.1:3]),'r');
hold off

%% 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
stairs(sort(Y1),[1:1000]/1000);
hold on
Y2=zeros(1,1000);
for k=1:1000
X=rpareto(5000,1.5);
Y2(k)=sqrt(5000)*(mean(X)-1/3);
end
stairs(sort(Y2),[1:1000]/1000,'r');
hold off
