import numpy as np
import pylab as plb

It=1000 #nombre d'iterations de l'algorithme de Metropolis Hastings pour
       

I=10 #nombre de types de sites
s=20 #population totale

lam=10*np.random.random(I) #1+10*rand(1,I,'u');//parametres des lois de Poisson

#initialisation des nombres de sites de chaque type occupes//

Y=np.floor(I*np.random.random(s)); #Tirage uniforme du type de site occupe par
                            #chacune des s particules
X=np.zeros(I) #vecteur des nombres de sites de chaque type occupes
for l in range(s):
  X[int(Y[l])]=X[int(Y[l])]+1 #comptage des nombres de sites de chaque type occupes



l=np.argmax(lam)  #indice de probabilite maximale dans la loi multinomiale
p=lam[l]/np.sum(lam) #probabilite correspondante

Z=np.zeros(It) #vecteur permettant de stocker X_k(l)

#Iterations de l'algorithme de Metropolis Hastings

for k in range(It):
  Z[k]=X[l]
  
  #####################################
  # ecrire le tirage d'un couple (i,j)#
  #####################################
  
  #################################################
  # avec probabilite d'acceptation donnee par le  #
  # ratio de Metropolis-Hastings, mettre a jour X #
  #################################################
 







# Pour le type de site l de probabilite maximale, trace en rouge n-->\sum_{k=0}^{n-1}X_k(l)/n
plb.plot(np.arange(It),np.cumsum(Z)/(1+np.arange(It)), color="r", label="frequence empirique")
# Trace en vert la limite theorique
plb.plot(np.arange(It),p*s*np.ones(It), color="g", label="limite")
# trace en bleu l'evolution de la moyenne empirique pour des tirages
# i.i.d. suivant la loi binomiale de X(l) sous la probabilite invariante
Y=1.*np.random.binomial(s,p,It)
plb.plot(np.arange(It),np.cumsum(Y)/(1+np.arange(It)), color="b", label="frequence empirique i.i.d.")

plb.show()