MHアルゴリズムの非対称提案分布のヘイスティング比g(x | x ')/ g(x' | x)を計算しますか?

Aug 18 2020

メトロポリスのアルゴリズムを理解しています。私が混乱するのは、非対称の提案分布が使用される可能性のあるMHアルゴリズムです。

P(x)とP(x ')は、ターゲット分布に従ったxとx'の尤度/確率密度を表すことを理解しています。同様に、g(x | x ')/ g(x' | x)は、非対称の提案分布を修正するために使用される用語であることを理解しています。私はその目的に混乱していません。その実行がわかりません。

トイプロブレムとして、指数分布サンプラーを開発しました。2つのバリエーションがあります。1つは対称的な提案分布を使用するもの、1つは一様分布です。そして、そうでないもの:すなわち、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

更新されたコードを見て、以下にプロットしてください。g(x | x ')/ g(x' | x)は、基本的に、として定義されている提案分布の下でこれらの摂動が見られる可能性の尺度であることに注意してくださいBeta(a=3,b=2) -0.5

まず、現在のイベントと提案されたイベントの違いを見つけます。次に、-0.5に調整します。これらの不偏摂動(-0.5はバイアス)と呼びます。3番目に各摂動の尤度を見つけます(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:]