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

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

//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);
kappa=zeros(M,1);//somme des valeurs de X(l) au cours des itérations

//boucle sur le nombre d'itération
n=it;
for nbit=1:4,
  for k=1:it,
    //tirage d'un couple (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);
    
    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;
    ///////////////////////
    //mettre à jour kappa//
    ///////////////////////
    
  end;



subplot(2,2,nbit)
/////////////////////////////////////////////////////////////////////////
//calculer moy=moyenne empirique de kappa/n et Z=moy+sqrt(n)(kappa/n-moy)
/////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////
//tracer l'histogramme de Z avec 20 bâtons//
////////////////////////////////////////////


v=variance(Z);
Z=sort(Z);
plot2d(Z,exp(-(Z-moy)^2/(2*v))/sqrt(2*%pi*v));
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;