%galton verifie la loi des grands nombres et le tcl
% pour une loi de reproduction choisie librement,
% et donne l'estimation de la probabilite d'extinction
% a partir de la fonction generatrice, ainsi qu'a partir
% d'une estimation empirique.
clf;
clear;
faute=1;
erreur='la somme n est pas egale a 1 une probabilite est negative';
while faute~=0,
p=input('Vecteur de probabilites ');
z=cumsum(p);
s=num2str(p);
faute=0;
if abs(z(length(z))-1)>eps faute=1;
end
if p~=abs(p) faute=2;
end
if faute~=0,
message=sprintf('Erreur! %s n est pas un vecteur de probabilites',s)
message=erreur(faute*30-29:faute*30)
end
end
pinv=p(length(p):-1:1);
pinvprime=(length(p)-1:-1:1).*pinv(1:length(p)-1);
v=p(length(p):-1:2);
uinv=cumsum(v);
moy=sum(uinv);
pext=pextinction(p,1,2);
n=input('Nombre de tirages de progenitures ');
x=rprogeniture(p,1,n);
for i=1:length(p),
b(i)=length(find(x==i-1))/n;
end
pextapp=pextinction(b,1,3);
a=zeros(length(p),n);
for i=1:n,
a(x(i)+1,i)=1;
end
figure(1);
l=lgn(a);
pause;
s=sprintf('Nombre d_echantillons de %d tirages de progenitures pour le tcl ',n);
m=input(s);
images=input('Nombre de pas de temps entre deux images consecutives ');
faute=1;
while faute==1,
pause(1);
s=sprintf('Quel nombre d_enfants voulez-vous etudier plus precisement entre 0 et %d ?',size(p)-1);
k=input(s);
kfinal=max(min(k,length(p)-1),0)+1;
if p(kfinal)*(1-p(kfinal))~=0
break;
end
message=sprintf('Erreur! on ne peut rien faire si p(%d)=%f',k,p(k))
end
faute=1;
while faute==1,
pas=input('Pas du dessin ');
xx0=input('Depart du dessin ');
xx1=input('Fin du dessin ');
faute=0;
if xx11,
plot(xxx,xyy,'r',xxx,xpp,'k');
else
plot(xxx,xpp,'k');
end
figure(1);
clf;
m=moviein(n);
for i=1:images:n,
clf;
xzz=pempirique(xtt(i,:),xxx);
plot(xxx,xyy,'r',xxx,xzz,'b');
m(:,i)=getframe;
end
%movie(m);
%extinction
clf;
clear
faute=1;
while faute==1,
p=input('Vecteur de probabilites ');
z=cumsum(p);
s=num2str(p);
faute=0;
if z(length(z))~=1 faute=1;
end
if p~=abs(p) faute=1;
end
if faute==1
message=sprintf('Erreur! %s n est pas un vecteur de probabilites',s)
end
pinv=p(length(p):-1:1);
end
pas=input('Pas du dessin de la fonction generatrice ');
iter=input('Nombre d iterations pour le dessin des fonctions generatrices ');
axis([0 1 0 1]);
x=0:pas:1;
y=x;
plot(y,x,'g');
hold on;
hold on;
x=polyval(pinv,x);
plot(y,x,'k');
hold on;
for i=1:iter-1,
x=polyval(pinv,x);
plot(y,x,'b');
hold on;
end
v=p(length(p):-1:2);
uinv=cumsum(v);
q=pinv;
q(length(q)-1)=q(length(q)-1)-1;
r=roots(q);
uu=cumsum(uinv);
m=uu(length(uu));
if m<=1,
message=sprintf('Comme l esperance de p vaut %f, la probabilite d extinction est egale a 1',m)
else
proba=r(find(r>0 & r<1));
message=sprintf('La probabilite d extinction vaut %f, et l esperance de p vaut %f',proba,m)
plot(proba,proba,'+r',proba,proba,'or');
hold on;
end
legende=sprintf('%d iterees successives de la fonction generatrice',iter-1);
legend('y=x','Fonction generatrice de la loi de reproduction',legende,2);
titre=sprintf('Determination de la probabilite d extinction pour la loi de reproduction p=%s',s);
title(titre);
function l=lgn(x)
%l=lgn(x) dessine l'evolution des moyennes successives des quantites
%associees a chaque ligne de la matrice x suivant le temps
%associe aux colonnes.
y=cumsum(x');
n=size(y,1);
v=1:n;
a=inv(diag(v,0));
l=a*y;
plot(l);
function d=pempirique(x,y)
%d=pempirique(x,y) retourne les valeurs de la fonction de
%repartition empirique de l'echantillon contenu dans x
%aux points de y. x et y sont des vecteurs unidimensionnels.
m=max(size(y));
n=max(size(x));
tx=sort(x);
for j=1:m,
s=max(find(y(j)>=tx));
if isempty(s)
s=0;
end
d(j)=s/n;
end
function proba=pextinction(p,dessins,ou);
%proba=pextinction(p,dessins,ou) determine la probablite d'extinction
% pour un processus de Galton-Watson de loi de reproduction
% p=[p_0 ... p_max], et dessine les courbes correspondantes
% lorsque dessins=1 dans la figure ou.
plettre=num2str(p);
pinv=p(length(p):-1:1);
if dessins==1
figure(ou);
clf;
pas=input('Pas du dessin de la fonction generatrice ');
iter=input('Nombre d iterations pour le dessin des fonctions generatrices ');
axis([0 1 0 1]);
x=0:pas:1;
y=x;
plot(y,x,'g');
hold on;
x=polyval(pinv,x);
plot(y,x,'k');
hold on;
for i=1:iter-1,
x=polyval(pinv,x);
plot(y,x,'b');
hold on;
end
legende=sprintf('%d iterees successives de la fonction generatrice',iter-1);
legend('y=x','Fonction generatrice de la loi de reproduction',legende,2);
titre=sprintf('Determination de la probabilite d extinction pour la loi de reproduction p=%s',plettre);
title(titre);
end
v=p(length(p):-1:2);
uinv=cumsum(v);
q=pinv;
q(length(q)-1)=q(length(q)-1)-1;
r=roots(q);
uu=cumsum(uinv);
moy=uu(length(uu));
if moy<=1,
proba=1;
message=sprintf('Comme l esperance de %s vaut %f, la probabilite d extinction est egale a 1',plettre,moy)
end
if moy>1,
probab=r(find(r>0 & r<1));
proba=min(probab);
message=sprintf('La probabilite d extinction vaut %f, et l esperance de p vaut %f',proba,moy)
end
if dessins==1,
plot(proba,proba,'+r',proba,proba,'or');
hold on;
end
function w=pprogeniture(p,x)
%w=pprogeniture(p,x) fournit une matrice p de meme taille
%que x des valeurs de la fonction de repartition de
%la loi p=(p_0,...,p_n) sur {0,...,n} en les points de x.
z=cumsum(p);
n=length(p);
a=size(x,1);
b=size(x,2);
for i=1:a,
for j=1:b,
k=min(floor(x(i,j)),n-1);
if k>=0 w(i,j)=z(k+1);
else w(i,j)=0;
end
end
end
function q=qprogeniture(p,x)
%q=qprogeniture(p,x) fournit une matrice q de meme taille
%que x des valeurs de la fonction de repartition inverse de
%la loi p=(p_0,...,p_n) sur {0,...,n} en les points de x.
z=cumsum(p);
n=length(p);
s=num2str(p);
a=size(x,1);
b=size(x,2);
faute=0;
message='';
if abs(z(n)-1)>eps faute=1;
end
if p~=abs(p) faute=1;
end
if faute==1
message=sprintf('Erreur! %s n est pas un vecteur de probabilites',s)
else
for i=1:a,
for j=1:b,
y=x(i,j);
if y<0
message=sprintf('Erreur! %f est negatif, je le remplace par 0',y)
y=0;
end
if y>1
message=sprintf('Erreur! %f est strictement superieur a 1, je le remplace par 1',y)
y=1;
end
k=1;
while y>z(k)
k=k+1;
if k>n
break;
end
end
q(i,j)=k-1;
end
end
end
function r=rprogeniture(p,a,b)
%r=rprogeniture(p,a,b) fournit une matrice r de taille (a,b) de nombres
%aleatoires suivant la loi p=(p_0,...,p_n)
if a<1
message=sprintf('%d est negatif ou nul!, je le remplace par 1',a)
a=1;
end
if b<1
message=sprintf('%d est negatif ou nul!, je le remplace par 1',b)
b=1;
end
x=rand(a,b);
r=qprogeniture(p,x);