////////////////////
//evenements rares//  
////////////////////

// Outil de calcul
function te=terme(k,l)
  te=exp(k*log(l)-sum(log([1:k])));
endfunction

// Proba pour qu'une poisson de parametre l depasse a
function pr=poissondep(a,l)
  k=ceil(a);
  te=terme(k,l);
  su=0;
  while (te>0);
    su=su+te;
    k=k+1;
    te=terme(k,l);
  end;
  pr=exp(-l)*su;
endfunction

// Comparaison des estimations de P(S_n>na) pour une somme de Poisson,
// bornes de Chebychev et de Chernov. Celle ci est optimale (3eme graphe)
function poissonrare(a,Nm,NM)
  ra=[Nm:NM]';
  ex=[];
  for n=Nm:NM;
    ex=[ex; poissondep(a*n,n)];
    //Alter: [P,Q]=cdfpoi("PQ",ceil(a*n)-1,n);
    //Alter: ex=[ex; Q];
  end;
  by=(a-1)^(-2)*ra.^(-1);
  Ia=1-a+a*log(a)
  rn=exp(-Ia*ra);
  //cv=((-log(ex)./ra)(NM-Nm+1)-Ia)
  xbasc();
  subplot(3,1,1);
  plot2d(ra,[ex by rn]);
  xtitle('Probabilites d''evenements rares pour une moyenne de variables de Poisson(1)');
  subplot(3,1,2);
  plot2d(ra,[ex rn], style=[1, 3]);
  subplot(3,1,3);
  plot2d(ra,[-log(ex)./ra Ia*ones(ra)], style=[1 5]);
  xtitle('La borne de Chernov est optimale a l''echelle logarithmique');
  legends(['valeurs theoriques' 'borne de Chebychev' 'borne de Chernov' ...
	   'valeur de la fonction de taux'],[1 2 3 5]);
endfunction

//////////////////////////
///// Galton Watson //////
//////////////////////////

//Inutile
//function al=tirage(pr)
//  u=rand();
//  P=cumsum(pr);
//  al=sum(1*(u>P));
//endfunction


// Simule une trajectoire de longueur n d'un Galton-Watson a generations
// binomiales(N,p)
// Attention, lorsque Np>1 (cas surcritique), ne pas prendre n trop
// grand, ca explose tres vite.
function galtwat(N,p,n)
  Z=[1];
  for i=1:n;
    Z=[Z sum(grand(1,Z(i),'bin',N,p))];
  end;
  xbasc();
  plot2d([0:n],Z,style=5);
  xtitle('Trajectoire d''un processus de Galton-Watson a generations binomiales ('+string(N)+','+ ...
	 string(p)+').');
endfunction

// Trajectoire a generations binomiale(N,p) stoppee en 5000. Rend 1 si
// extinction, 0 sinon.
function ext=galtwat2(N,p)
  Z=[1];
  i=1;
  while ((Z(i)>0) & (Z(i)<5000));
    Z=[Z sum(grand(1,Z(i),'bin',N,p))];
    i=i+1;
  end;
  ext=1*(Z(i)==0);
  //xbasc();
  //plot2d([1:length(Z)],Z);
endfunction

// Calcul (approche) de pi et G'(pi) pour des generations binomiales (N,p)
function [cal,der]=calpi(N,p)
  cp=1;
  cs=0;
  while (abs(cp-cs)>0.0000000001);
    cp=cs;
    cs=((1-p)+p*cp)^N;
  end;
  cal=cs;
  der=N*p*(1-p+p*cs)^(N-1);
endfunction

// Estimation de pi par Monte Carlo dans le cas binomial. Comparaison
// avec la valeur calculee
// Rend valeur estimee finale et valeur calculee
function [est,cal]=estpi(N,p,m)
  val=[];
  for i=1:m;
    val=[val;galtwat2(N,p)];
  end;
  esti=cumsum(val)./[1:m]';
  [cal,der]=calpi(N,p);
  est=esti(m);
  xbasc();
  plot2d([1:m]',[esti esti-1.96*sqrt(esti.*(1-esti))./sqrt([1:m]') ...
		 esti+1.96*sqrt(esti.*(1-esti))./sqrt([1:m]') cal* ...
		 ones([1:m]')],style=[1,2,2,5]);
  xtitle('Estimation de pi pour une loi binomiale ('+string(N)+','+ ...
	 string(p)+').');
  legends(['Valeur calculee','Valeur estimee','Intervalle de confiance'],[5,1,2]);
endfunction

// Comme galtwat2 pour une geometrique
function ext=galtwatgeom(p)
  Z=[1];
  i=1;
  while ((Z(i)>0) & (Z(i)<5000));
    Z=[Z sum(grand(1,Z(i),'nbn',1,p))];
    i=i+1;
  end;
  ext=1*(Z(i)==0);
  //xbasc();
  //plot2d([1:length(Z)],Z);
endfunction
 
// Comme estpi pour une geometrique. On connait alors explicitement pi
function est=estpigeom(p,m)
  val=[];
  for i=1:m;
    val=[val;galtwatgeom(p)];
  end;
  esti=cumsum(val)./[1:m]';
  est=esti(m);
  xbasc();
  l=min(1,p/(1-p));
  plot2d([1:m]',[esti esti-1.96*sqrt(esti.*(1-esti))./sqrt([1:m]') ...
		 esti+1.96*sqrt(esti.*(1-esti))./sqrt([1:m]') l*ones([1:m]')],style=[1,2,2,5]);
endfunction

// Variante de galtwat2 qui s'arrete a l'etape n ou au seuil a+5000.
// Rend la trajectoire complete (et meme completee a a+5000 si ca depasse
function Z=galtwat3(N,p,n,a)
  Z=[1];
  i=1;
  while (i<=n & (Z(i)<(a+5000)));
    Z=[Z sum(grand(1,Z(i),'bin',N,p))];
    i=i+1;
  end;
  if (length(Z)<(n+1)) then
    Z=[Z (a+5000)*ones(1,n+1-length(Z))];
  end;
  endfunction

// Illustre l'explosion exponentielle dans le cas surcritique. On trace 
// log(1-P(Z_n>A<|non extinction)) et on compare a la droite de pente 
// log(G'(pi))
function surcritique(N,p,m,n,A)
  est=[];
  for i=1:length(A);
    esta=[];
    a=A(i);
    Z=[];
    k=0;
    for j=1:m;
      z=galtwat3(N,p,n,a);
      if (z(n+1)>0) then 
	Z=[Z;z]; 
	k=k+1;
      end;
    end;
    est=[est log(1-(sum(1*(Z>a),1))'/k)];
  end;
  [cal,der]=calpi(N,p);
  xbasc();
  plot2d([0:n]',[log(der)*[0:n]' est],style=[5 1:3]);
  xtitle('Explosion exponentielle dans le cas surcritique (echelle logarithmique)');
  legends(['Pente log(G''(pi))','a='+string(A(1))+'.','a='+string(A(2))+ ...
	   '.','a='+string(A(3))+'.'],[ 5 1:3]);
endfunction