La méthode MCMC : Markoc Chain Monte Carlo (cours)

[@menadhara1352]

bonjours,
question: comment utiliser q pour généré x ?
merci pour votre réponse.

[@Achillionable]

Un moyen est en fait d'utiliser une sorte de Monte-Carlo pour ça.
C'est brute force, mais ça marche
Un petit code python pour illustrer :
import numpy as np
import matplotlib.pyplot as plt
X=np.linspace(0,5,1000)
def Y(X): ## ma distribution q pour obtenir x
return(np.sin(X*np.pi/5)*np.abs(np.sin(1/3*X**2)))
from scipy.integrate import quad
I=quad(Y,0,5,args=(X)) ## j'en fait l'intégrale pour être sûr d'avoir une pdf avec intégrale=1
def Y(X): ## ma distribution q pour obtenir x, intégrale=1
return(np.sin(X*np.pi/5)/I*np.abs(np.sin(1/3*X**2)))
L=[]
for _ in range(10000000):
x=np.random.rand()*5
if Y(x)>np.random.rand(): ## L'acceptation, pour avoir le nouveau x à partir de q
L+=[x]
H=np.histogram(L,bins=50)
plt.hist(L,bins=50)
plt.plot(X,Y(X)/np.max(Y(X))*np.max(H[0])) ## Je vérifie que mes tirages suivent en effet la distribution q.
Et on observe que ça marche bien en effet.
Du coup, à partir d'un q classique (chez moi Y), il suffit de faire ce petit tirage. ça fonctionnera forcément car la distribution q (ici Y) est une pdf, donc elle ne peut être supérieure à 1. Le plus efficace serait de vérifier que Y(x)> a ϵ [0,max(Y)], pour éviter les tirages inutiles.
maxi=np.max(Y(X))
x=np.random.rand()*5
if Y(x)>np.random.rand()*maxi
Et le prochain x' sera l'x qui respecte la condition if.
Il y a probablement plus simple, mais ça fonctionne, j'ai directement pensé à ça.