n=10;
N=1000;
lambda=3/4;
mu=[1/3 1/2 1/6];
xbasc();

//simulation suivant la loi de Poisson de paramètre lambda
function [z]=poissonruben(v);
  k=-1;
  b=0;
  t=0;
  while (b==0)
    u=rand();
    t=t-log(u)/lambda;
    if (t>1)
      b=1;
    end,
    k=k+1;
  end,
  z=k;
endfunction

//simulation suivant la loi d'un vecteur ligne
function [z]=musim(u);
  //mu doit être un vecteur ligne
 k=length(mu);
 P=mu(1:k-1);
 x=ones(1,k);
 y=cumsum(x);
 z=y*grand(1,'mul',1,P')-1;
endfunction

//simulation d'un accroissement de Galton-Watson
function [z]=acc1(x)
  s=0;
  for i=1:x
    s=s+poisson(1);
  end,
  z=s;
endfunction

//simulation d'un accroissement de Galton-Watson
function [z]=acc2(x)
  s=0;
  for i=1:x
    s=s+musim(1);
  end,
  z=s;
endfunction


X1=ones(1,N);
X2=ones(1,N);
z1=zeros(1,n+1);
z2=zeros(1,n+1);
l1=zeros(1,n+1);
for i=1:n
  Y1=feval(X1,acc1);
  Y2=feval(X2,acc2);
  X1=Y1;
  X2=Y2;
  w1=find(X1>0);
  w2=find(X2>0);
  k1=length(w1);
  k2=length(w2);
  l1(1,i+1)=k1;
  z1(1,i+1)=log(k1/N);
  z2(1,i+1)=log(k2/N);
end,  
t=[0:1:n];

xset("windows",0);
plot(t,z1,'red');
plot(t,log(3/4)*t,'yellow');
xset("windows",1);
plot(t,z2,'blue');
plot(t,log(0.83)*t,'cyan');



 
