Calcul du rapport de hastings, g (x | x ') / g (x' | x) pour les distributions de proposition asymétriques dans l'algorithme MH?
Je comprends l'algorithme Metropolis. Là où je suis confus, c'est l'algorithme MH où des distributions de proposition asymétriques peuvent être utilisées.
Je comprends que P (x) et P (x ') représentent la densité de vraisemblance / probabilité de x et x' selon la distribution cible. De même, je comprends que g (x | x ') / g (x' | x) est un terme utilisé pour corriger une distribution de proposition asymétrique. Je ne suis pas confus par son objectif; Je ne comprends pas son exécution.
En tant que problème de jouet, j'ai développé un échantillonneur de distribution exponentielle. Il existe deux variantes, l'une qui utilise une distribution de proposition symétrique, l'uniforme dist. Et celui qui ne le fait pas: à savoir, Beta(a=3,b=2) - 0.5
. J'ai choisi cette distribution parce que (A) elle est asymétrique et principalement positive (mais parfois négative, en raison du terme -0,5.)
Je ne sais pas comment trouver g(x|x')/g(x'|x)
.
Code:
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = lam
for i in range(hops):
states.append(current)
# movement = current + random.uniform(-1,1) # does not require asymmetric correction
movement = current + np.random.beta(a=3,b=2)-0.5 # requires asymmetric correction
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
lam = 1
exp_samples = exponential_MH(hops=10_000,lam=lam)
lines = np.linspace(0,5,10_000)
exp_curve = [lam*np.exp(-l*lam) for l in lines]
plt.hist(exp_samples,normed=1,bins=20)
plt.plot(lines,exp_curve)
Diagramme de distribution de proposition asymétrique:
Diagramme de distribution de proposition symétrique (avec les lignes commentées inversées):
Pour répondre à cette question, veuillez modifier le code afin de refléter une g(x|x')/g(x'|x)
distribution valide pour la proposition Beta(a=3,b=2) - 0.5
dont les perturbations sont tirées.
Réponses
Jetez un œil au code mis à jour et au graphique ci-dessous. Notez que g (x | x ') / g (x' | x) est essentiellement une mesure de la probabilité que ces perturbations soient vues sous la distribution de proposition, qui a été définie comme Beta(a=3,b=2) -0.5
.
Tout d'abord, trouvez la différence entre les événements actuels et proposés. Deuxièmement, ajustez pour le -0,5; nous appellerons ces perturbations non biaisées (où -0,5 est un biais.) en troisième lieu, trouver la probabilité de chaque perturbation (curr-> prop & prop-> curr). Enfin, renvoyez le ratio sous la forme correction
.
Nous allons utiliser cette correction et la multiplier par d'autres termes dans la acceptance
définition de variable. C'est à peu près tout!
def target(x,lam):
return int(x>0) * lam * np.exp(-x * lam)
def correct(prop,curr,a=3,b=2):
x0 = curr - prop + 0.5
x1 = prop - curr + 0.5
b0 = beta.pdf(x=x0, a=a, b=b)
b1 = beta.pdf(x=x1, a=a, b=b)
return b0/b1
def exponential_MH(hops,lam=3):
states = []
burn_in = int(hops*0.2)
current = 1
for i in range(hops):
states.append(current)
movement = current + np.random.beta(a=3,b=2)-0.5 # requires assymetric correction
correction = correct(curr=current,prop=movement)
acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)*correction
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]