Menghitung rasio hastings, g (x | x ') / g (x' | x) untuk distribusi proposal asimetris dalam algoritma MH?
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.5
yang berasal dari gangguan.
Jawaban
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 acceptance
definisi 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:]
