MH algoritmasında asimetrik öneri dağılımları için hastings oranı, g (x | x ') / g (x' | x) hesaplanıyor mu?

Question

Metropolis algoritmasını anlıyorum. Kafamın karıştığı yer, asimetrik teklif dağılımlarının kullanılabileceği MH algoritmasıdır.

Hedef dağılıma göre P (x) ve P (x ') x ve x' olasılık / olasılık yoğunluğunu temsil ettiğini anlıyorum. Aynı şekilde, g (x | x ') / g (x' | x) 'in asimetrik bir teklif dağılımını düzeltmek için kullanılan bir terim olduğunu anlıyorum. Amacı kafam karışmıyor; İnfazını anlamıyorum.

Bir oyuncak problemi olarak üstel dağılım örnekleyicisi geliştirdim. Simetrik bir öneri dağılımını kullanan iki varyasyon vardır, tekbiçimli uzaklık. Ve yapmayan: Yani Beta(a=3,b=2) - 0.5,. Bu dağılımı seçtim çünkü (A) asimetrik ve çoğunlukla pozitif (ancak -0.5 terimi nedeniyle bazen negatif).

Nasıl bulacağımı bilmiyorum g(x|x')/g(x'|x).

Kod:

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)

Asimetrik teklif dağıtım grafiği:

Simetrik teklif dağıtım grafiği (yorumlu satırlar çevrilmiş olarak):

Bu soruyu cevaplamak için, lütfen kodu, tedirginliklerin g(x|x')/g(x'|x)alındığı teklif dağıtımı için geçerli olacak şekilde düzenleyin Beta(a=3,b=2) - 0.5.

Answer 1

Aşağıdaki güncellenmiş koda ve arsaya bir göz atın. G (x | x ') / g (x' | x) 'in esasen, olarak tanımlanan teklif dağıtımı altında bu tedirginliklerin ne kadar olası olduğunun bir ölçüsü olduğuna dikkat edin Beta(a=3,b=2) -0.5.

İlk olarak, mevcut ve önerilen olaylar arasındaki farkı bulun. İkinci olarak, -0,5 için ayarlayın; bu tarafsız tedirginlikler adını vereceğiz (burada -0.5 bir önyargıdır.) üçüncü olarak her bir pertürbasyonun olasılığını bulun (akım-> prop & prop-> akım). Son olarak, oranı olarak döndürün correction.

Bu düzeltmeyi kullanacağız ve acceptancedeğişken tanımındaki diğer terimlerle çarpacağız. Hepsi bukadar!

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