Calcul du rapport de hastings, g (x | x ') / g (x' | x) pour les distributions de proposition asymétriques dans l'algorithme MH?

Aug 18 2020

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.5dont les perturbations sont tirées.

Réponses

2 jbuddy_13 Aug 18 2020 at 02:19

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 acceptancedé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:]