एमएच एल्गोरिथ्म में असममित प्रस्ताव वितरण के लिए हेस्टिंग्स अनुपात, जी (x | x ') / g (x' | x) की गणना करना।
मैं मेट्रोपोलिस एल्गोरिथ्म को समझता हूं। जहां मैं भ्रमित हो जाता हूं वह एमएच एल्गोरिथम है जहां असममित प्रस्ताव वितरण का उपयोग किया जा सकता है।
मैं समझता हूं कि पी (एक्स) और पी (एक्स ') लक्ष्य वितरण के अनुसार एक्स और एक्स' की संभावना / संभावना घनत्व का प्रतिनिधित्व करते हैं। इसी तरह, मैं समझता हूँ कि g (x | x ') / g (x' | x) एक शब्द है जिसका उपयोग विषम प्रस्ताव वितरण को सही करने के लिए किया जाता है। मैं इसके उद्देश्य से भ्रमित नहीं हूँ; मैं इसके निष्पादन को नहीं समझता।
एक खिलौना समस्या के रूप में, मैंने एक घातीय वितरण नमूना विकसित किया है। दो भिन्नताएं हैं, एक जो सममित प्रस्ताव वितरण का उपयोग करती है, एक समान डिस्टर्ब। और एक है कि नहीं: अर्थात् Beta(a=3,b=2) - 0.5। मैंने इस वितरण को इसलिए चुना है क्योंकि (ए) यह असममित है और ज्यादातर सकारात्मक है (लेकिन कभी-कभी नकारात्मक, -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जिसे वितरण से खींचा गया है।
जवाब
नीचे दिए गए अद्यतन कोड और प्लॉट पर एक नज़र डालें। ध्यान दें कि g (x | x ') / g (x' | x) अनिवार्य रूप से इस बात का एक उपाय है कि प्रस्ताव वितरण के तहत इन गड़बड़ियों की कितनी संभावना है, जिसे परिभाषित किया गया है Beta(a=3,b=2) -0.5।
सबसे पहले, वर्तमान और प्रस्तावित घटनाओं के बीच अंतर खोजें। दूसरा, -0.5 के लिए समायोजित करें; हम इन निष्पक्ष गड़बड़ियों (जहां -0.5 एक पूर्वाग्रह है) को बुलाएंगे। तीसरा प्रत्येक गड़बड़ी (वक्र-> प्रोप और प्रोप-> वक्र) की संभावना का पता लगाएं। अंत में, इस अनुपात को वापस लौटाएं 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:]
