clf;
M=10000//taille échantillon Monte Carlo
I=10;//nombre de types de sites
s=20;//population totale
it=10;//nombre initial d'itérations


lam=10*rand(1,I,'u');//1+10*rand(1,I,'u');//paramètres des lois de Poisson
Lam=sum(lam);

//initialisation des nombres de sites de chaque type occupés//


Y=1+floor(I*rand(M,s,'u')); //Tirage uniforme du type de site occupé par
                            //chacune des s particules
X=zeros(M,I);//vecteur des nombres de sites de chaque type occupés
for m=1:M,
  for l=1:s,
    X(m,Y(m,l))=X(m,Y(m,l))+1; //comptage des nombres de sites de chaque type occupés
  end;
end;
//X=[20*ones(M,1),zeros(M,I-1)];//initialisation déterministe

[ma,l]=max(lam);

n=it
//boucle sur le nombre d'itérations
for nbit=1:4,
  for k=1:it,
    //tirage de M couples (i,j)
    i=1+sum(cumsum(X,'c')<s*rand(M,1)*ones(1:I),'c');
    j=1+floor((I-1)*rand(M,1));
    j=j+(j>=i);
    //mise à jour des M états 
    inc=(rand(1,M)<lam(j)./lam(i))';
    i=(1:M)'+M*(i-1);    
    X(i)=X(i)-inc;
    j=(1:M)'+M*(j-1);
    X(j)=X(j)+inc;
  end;

  compt=zeros(0:s);
  //comptage dans compt(x+1) du nombre de simulations m pour lesquelles
  //X(m,l)=x
  for m=1:M,
    compt(X(m,l)+1)=compt(X(m,l)+1)+1; 
  end;
  compt/M

  
  subplot(2,2,nbit)// Trace le diagramme en bâton des moyennes empiriques des
  // proportions de particules dans les types de sites
  plot2d3((0:s), compt/M,rect=[-1,0,s+1,max(compt/M)*1.1],nax=[0,s+3,1,9]);
  // Trace en bleu la loi binomiale
  plot2d3((0:s)+0.2,binomial(lam(l)/Lam,s),style=4); 
  xtitle('n='+string(n))
it=n*3;//multiplication par 3 du nombre d'itérations supplémentaires
n=n*4; //multiplication par 4 de l'horizon 
end;