clf;
It=1000;//nombre d'itérations de l'algorithme de Métropolis Hastings pour
       

I=10;//nombre de types de sites
s=20;//population totale

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(1,s,'u')); //Tirage uniforme du type de site occupé par
                            //chacune des s particules
X=zeros(1:I);//vecteur des nombres de sites de chaque type occupés
for l=1:s,
  X(Y(l))=X(Y(l))+1; //comptage des nombres de sites de chaque type occupés
end;


[m,l]=max(lam);//indice de probabilité maximale dans la loi multinomiale
Z=[];//vecteur permettant de stocker X_k(l)

//Itérations de l'algorithme de Métropolis Hastings

for k=1:It,
  //tirage d'un couple (i,j)//
  i=1+sum(cumsum(X)<s*rand(1));
  j=1+floor((I-1)*rand(1));
  if (j>=i)
    j=j+1;
  end;
  if (rand(1)<lam(j)/lam(i))//cas Poissonien
    //Cas binomial (rand(1)<p(j)*(1-p(i))*(m(j)-X(j))/(p(i)*(1-p(j))*(m(i)+1-X(i))))
    X(i)=X(i)-1;
    X(j)=X(j)+1;
  end;
  Z=[Z,X(l)];
end;





p=lam(l)/sum(lam);


// Pour le type de site l de probabilité maximale, trace en n-->\sum_{k=0}^{n-1}X_k(l)/n
plot2d((1:It),cumsum(Z)./(1:It),style=5)
// Trace en vert la limite théorique
plot2d((1:It),20*p*ones(1:It), style=3);
// trace en bleu l'évolution de la moyenne empirique pour des tirages
// i.i.d. suivant la loi binomiale de X(l) sous la probabilité invariante
Y=grand(1,It,'bin',s,p);
plot2d((1:It),cumsum(Y)./(1:It),style=2)