Menghitung rasio hastings, g (x | x ') / g (x' | x) untuk distribusi proposal asimetris dalam algoritma MH?

Aug 18 2020

Saya memahami algoritma Metropolis. Di mana saya bingung adalah algoritma MH di mana distribusi proposal asimetris dapat digunakan.

Saya memahami bahwa P (x) dan P (x ') mewakili kepadatan kemungkinan / probabilitas dari x dan x' sesuai dengan distribusi target. Demikian pula, saya memahami bahwa g (x | x ') / g (x' | x) adalah istilah yang digunakan untuk mengoreksi distribusi proposal asimetris. Saya tidak bingung dengan tujuannya; Saya tidak mengerti pelaksanaannya.

Sebagai masalah mainan, saya telah mengembangkan sampler distribusi eksponensial. Ada dua variasi, yang menggunakan distribusi proposal simetris, dist seragam. Dan satu yang tidak: Yaitu , Beta(a=3,b=2) - 0.5. Saya memilih distribusi ini karena (A) asimetris dan sebagian besar positif (tetapi terkadang negatif, karena suku -0,5.)

Saya tidak tahu bagaimana menemukannya g(x|x')/g(x'|x).

Kode:

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)

Plot distribusi proposal asimetris:

Plot distribusi proposal simetris (dengan baris komentar dibalik):

Untuk menjawab pertanyaan ini, harap edit kode untuk mencerminkan validitas g(x|x')/g(x'|x)distribusi proposal Beta(a=3,b=2) - 0.5yang berasal dari gangguan.

Jawaban

2 jbuddy_13 Aug 18 2020 at 02:19

Lihatlah kode dan plot yang diperbarui di bawah ini. Perhatikan bahwa g (x | x ') / g (x' | x) pada dasarnya adalah ukuran seberapa besar kemungkinan gangguan ini terlihat dalam distribusi proposal, yang telah didefinisikan sebagai Beta(a=3,b=2) -0.5.

Pertama, temukan perbedaan antara peristiwa saat ini dan yang diusulkan. Kedua, sesuaikan untuk -0,5; kita akan menyebut gangguan yang tidak bias ini (di mana -0.5 adalah bias.) ketiga temukan kemungkinan dari setiap gangguan (arus-> prop & prop-> arus). Terakhir, kembalikan rasio sebagai correction.

Kami akan menggunakan koreksi ini dan mengalikannya dengan istilah lain dalam acceptancedefinisi variabel. Cukup banyak!

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