MH 알고리즘에서 비대칭 제안 분포에 대한 헤이스팅 비율 g (x | x ') / g (x'| x) 계산?

Question

Metropolis 알고리즘을 이해합니다. 내가 헷갈리는 부분은 비대칭 제안 분포가 사용될 수있는 MH 알고리즘입니다.

나는 P (x)와 P (x ')가 목표 분포에 따라 x와 x'의 우도 / 확률 밀도를 나타낸다는 것을 이해합니다. 마찬가지로 g (x | x ') / g (x'| x)가 비대칭 제안 분포를 수정하는 데 사용되는 용어라는 것을 이해합니다. 나는 그것의 목적으로 혼동하지 않는다. 나는 그것의 실행을 이해하지 못한다.

장난감 문제로 지수 분포 샘플러를 개발했습니다. 두 가지 변형이 있습니다. 하나는 대칭 제안 분포 인 uniform dist를 사용합니다. 그리고 그렇지 않은 사람 : 즉, 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.

Answer 1

아래의 업데이트 된 코드와 플롯을 살펴보십시오. g (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:]