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

Aug 18 2020

ฉันเข้าใจอัลกอริทึม 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ที่มาจากการรบกวน

คำตอบ

2 jbuddy_13 Aug 18 2020 at 02:19

ดูโค้ดที่อัปเดตและพล็อตด้านล่าง ขอให้สังเกตว่ากรัม (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:]