clf;
I=10;//nombre de types de sites
s=20;//population totale
It=1000;//nombre d'itérations de l'algorithme de Métropolis Hastings pour
       



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,
  //////////////////////////////////////
  //écrire le tirage d'un couple (i,j)//
  //////////////////////////////////////
  
  
  
  ////////////////////////////////////////////////////////////////////
  //écrire la mise à jour de X avec la probabilité d'acceptation donnée
  //par le rapport de Metropolis Hastings
  /////////////////////////////////////////////////////////////////////
  
  
  
  
  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))
// 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);