การคำนวณอัตราส่วน hastings, g (x | x ') / g (x' | x) สำหรับการแจกแจงข้อเสนอแบบไม่สมมาตรในอัลกอริทึม MH?
ฉันเข้าใจอัลกอริทึม Metropolis ที่ฉันสับสนคืออัลกอริทึม MH ที่อาจใช้การแจกแจงข้อเสนอแบบไม่สมมาตร

ฉันเข้าใจว่า P (x) และ P (x ') แสดงถึงความเป็นไปได้ / ความหนาแน่นของความน่าจะเป็นของ x และ x' ตามการแจกแจงเป้าหมาย ในทำนองเดียวกันฉันเข้าใจว่า g (x | x ') / g (x' | x) เป็นคำที่ใช้เพื่อแก้ไขการแจกแจงข้อเสนอแบบอสมมาตร ฉันไม่สับสนกับจุดประสงค์ของมัน ฉันไม่เข้าใจการทำงานของมัน
ในฐานะที่เป็นปัญหาเกี่ยวกับของเล่นฉันได้พัฒนาตัวสุ่มตัวอย่างการกระจายเลขชี้กำลัง มีสองรูปแบบหนึ่งที่ใช้การแจกแจงข้อเสนอแบบสมมาตรระยะสม่ำเสมอ และสิ่งที่ไม่: ได้แก่Beta(a=3,b=2) - 0.5
. ฉันเลือกการแจกแจงแบบนี้เพราะ (A) มันไม่สมมาตรและส่วนใหญ่เป็นค่าบวก (แต่บางครั้งก็เป็นลบเนื่องจากเทอม -0.5)
ฉันไม่รู้ว่าจะหาได้g(x|x')/g(x'|x)
อย่างไร
รหัส:
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)
แผนผังการกระจายข้อเสนอแบบไม่สมมาตร:

โครงร่างการกระจายข้อเสนอแบบสมมาตร (โดยพลิกบรรทัดที่แสดงความคิดเห็น):

หากต้องการตอบคำถามนี้โปรดแก้ไขโค้ดเพื่อให้แสดงถึงความถูกต้องg(x|x')/g(x'|x)
สำหรับการแจกจ่ายข้อเสนอBeta(a=3,b=2) - 0.5
ที่มาจากการรบกวน
คำตอบ
ดูโค้ดที่อัปเดตและพล็อตด้านล่าง ขอให้สังเกตว่ากรัม (x | x ') / g (x' | x) Beta(a=3,b=2) -0.5
เป็นหลักวัดว่าแนวโน้มเยี่ยงอย่างเหล่านี้จะได้รับการเห็นภายใต้ข้อเสนอการจัดจำหน่ายที่ได้รับการกำหนดให้เป็น
ขั้นแรกค้นหาความแตกต่างระหว่างเหตุการณ์ปัจจุบันและเหตุการณ์ที่เสนอ ประการที่สองปรับสำหรับ -0.5; เราจะเรียกการก่อกวนที่ไม่เอนเอียงเหล่านี้ (โดยที่ -0.5 คืออคติ) ประการที่สามค้นหาความเป็นไปได้ของการก่อกวนแต่ละครั้ง (curr-> prop & prop-> curr) สุดท้ายคืนอัตราส่วนเป็นcorrection
.
เราจะใช้การแก้ไขนี้และคูณกับคำอื่น ๆ ในacceptance
นิยามตัวแปร มันสวยมาก!
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:]
