Calcolando il rapporto hastings, g (x | x ') / g (x' | x) per distribuzioni di proposte asimmetriche nell'algoritmo MH?

Aug 18 2020

Capisco l'algoritmo di Metropolis. Dove mi confondo è l'algoritmo MH in cui possono essere utilizzate distribuzioni di proposte asimmetriche.

Capisco che P (x) e P (x ') rappresentano la probabilità / densità di probabilità di x e x' in base alla distribuzione target. Allo stesso modo, capisco che g (x | x ') / g (x' | x) è un termine usato per correggere una distribuzione asimmetrica della proposta. Non sono confuso dal suo scopo; Non capisco la sua esecuzione.

Come problema di giocattoli, ho sviluppato un campionatore di distribuzione esponenziale. Esistono due varianti, una che utilizza una distribuzione proposta simmetrica, la dist. Uniforme. E uno che non lo fa: Vale a dire, Beta(a=3,b=2) - 0.5. Ho scelto questa distribuzione perché (A) è asimmetrica e per lo più positiva (ma occasionalmente negativa, a causa del termine -0,5).

Non ho idea di come trovarlo g(x|x')/g(x'|x).

Codice:

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)

Grafico di distribuzione delle proposte asimmetriche:

Grafico di distribuzione simmetrico della proposta (con righe commentate invertite):

Per rispondere a questa domanda, modifica il codice in modo che rifletta un valido g(x|x')/g(x'|x)per la distribuzione della proposta da Beta(a=3,b=2) - 0.5cui derivano le perturbazioni.

Risposte

2 jbuddy_13 Aug 18 2020 at 02:19

Dai un'occhiata al codice aggiornato e alla trama di seguito. Si noti che g (x | x ') / g (x' | x) è essenzialmente una misura della probabilità che queste perturbazioni siano viste sotto la distribuzione della proposta, che è stata definita come Beta(a=3,b=2) -0.5.

Innanzitutto, trova la differenza tra gli eventi attuali e quelli proposti. Secondo, regolare per il -0,5; chiameremo queste perturbazioni non distorte (dove -0.5 è un bias.) terzo trovare la probabilità di ogni perturbazione (curr-> prop & prop-> curr). Infine, restituisci il rapporto come correction.

Useremo questa correzione e la moltiplicheremo con altri termini nella acceptancedefinizione della variabile. Questo è praticamente tutto!

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