了解具有不對稱提案分佈的 Metropolis-Hastings

July 14, 2013

我一直在嘗試理解 Metropolis-Hastings 算法，以便編寫代碼來估計模型的參數（即）。根據參考書目，Metropolis-Hastings 算法具有以下步驟：

產生

在哪裡

我想問幾個問題：

參考書目指出，如果是一個對稱分佈的比率變為 1，算法稱為 Metropolis。那是對的嗎？Metropolis 和 Metropolis-Hastings 之間的唯一區別是第一個使用對稱分佈？“隨機漫步”大都會（-黑斯廷斯）怎麼樣？它與其他兩個有何不同？

我在網上找到的大多數示例代碼都使用高斯提議分佈因此在哪裡是似然比。如果提案分佈是泊松分佈怎麼辦？我想我理性地理解為什麼在使用非對稱分佈時該比率不會變為 1，但我不太確定是否在數學上理解它或如何用代碼實現它。有人可以給我一個使用非對稱提議分佈的 Metropolis-Hastings 算法的簡單代碼（C、python、R、偽代碼或任何你喜歡的）示例嗎？

參考書目指出，如果 q 是對稱分佈，則比率 q(x|y)/q(y|x) 變為 1，該算法稱為 Metropolis。那是對的嗎？

是的，這是正確的。Metropolis 算法是 MH 算法的一個特例。

“隨機漫步”大都會（-黑斯廷斯）怎麼樣？它與其他兩個有何不同？

在隨機遊走中，提案分佈在每一步之後重新以鏈最後生成的值為中心。通常，在隨機遊走中，提議分佈是高斯分佈，在這種情況下，這種隨機遊走滿足對稱性要求，算法是 Metropolis。我想您可以執行具有不對稱分佈的“偽”隨機遊走，這會導致提案在偏斜的相反方向上過於漂移（左偏分佈會偏向右側的提案）。我不知道你為什麼要這樣做，但你可以而且這將是一個大都會黑斯廷斯算法（即需要額外的比率項）。

它與其他兩個有何不同？

在非隨機遊走算法中，提議分佈是固定的。在隨機遊走變體中，提案分佈的中心在每次迭代時都會發生變化。

如果提案分佈是泊松分佈怎麼辦？

然後你需要使用 MH 而不僅僅是 Metropolis。大概這將是對離散分佈進行採樣，否則您將不想使用離散函數來生成您的建議。

在任何情況下，如果採樣分佈被截斷或者您事先知道它的偏斜，您可能想要使用非對稱採樣分佈，因此需要使用 Metropolis-hastings。

有人可以給我一個簡單的代碼（C、python、R、偽代碼或任何你喜歡的）示例嗎？

這裡是大都會：
Metropolis <- function(F_sample # distribution we want to sample
                     , F_prop  # proposal distribution 
                     , I=1e5   # iterations
              ){
 y = rep(NA,T)
 y[1] = 0    # starting location for random walk
 accepted = c(1)

 for(t in 2:I)    {
   #y.prop <- rnorm(1, y[t-1], sqrt(sigma2) ) # random walk proposal
   y.prop <- F_prop(y[t-1]) # implementation assumes a random walk. 
                            # discard this input for a fixed proposal distribution

   # We work with the log-likelihoods for numeric stability.
   logR = sum(log(F_sample(y.prop))) -
          sum(log(F_sample(y[t-1])))    

   R = exp(logR)

   u <- runif(1)        ## uniform variable to determine acceptance
   if(u < R){           ## accept the new value
     y[t] = y.prop
     accepted = c(accepted,1)
   }    
   else{
     y[t] = y[t-1]      ## reject the new value
     accepted = c(accepted,0)
   }    
 }
 return(list(y, accepted))
}
讓我們嘗試使用它來採樣雙峰分佈。首先，讓我們看看如果我們對 propsal 使用隨機遊走會發生什麼：
set.seed(100)

test = function(x){dnorm(x,-5,1)+dnorm(x,7,3)}

# random walk
response1 <- Metropolis(F_sample = test
                      ,F_prop = function(x){rnorm(1, x, sqrt(0.5) )}
                     ,I=1e5
                      )
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
# 0.85585   not bad

# looks about how we'd expect
plot(density(y_trace1))
abline(v=-5);abline(v=7) # Highlight the approximate modes of the true distribution
現在讓我們嘗試使用固定的提案分佈進行抽樣，看看會發生什麼：
response2 <- Metropolis(F_sample = test
                           ,F_prop = function(x){rnorm(1, -5, sqrt(0.5) )}
                           ,I=1e5
                      )

y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # .871, not bad
起初這看起來不錯，但如果我們看一下後密度……
plot(density(y_trace2))
我們會看到它完全停留在局部最大值。這並不完全令人驚訝，因為我們實際上將我們的提案分佈集中在那裡。如果我們將其集中在另一種模式上，也會發生同樣的事情：
response2b <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, 7, sqrt(10) )}
                       ,I=1e5
)

plot(density(response2b[[1]]))
我們可以嘗試在兩種模式之間放棄我們的建議，但我們需要將方差設置得非常高，以便有機會探索其中的任何一種
response3 <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, -2, sqrt(10) )}
                       ,I=1e5
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # .3958! 
請注意，我們的提案分佈中心的選擇如何對我們的採樣器的接受率產生重大影響。
plot(density(y_trace3))
plot(y_trace3) # we really need to set the variance pretty high to catch 
              # the mode at +7. We're still just barely exploring it
我們仍然陷入兩種模式中的較近者。讓我們嘗試在兩種模式之間直接刪除它。
response4 <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, 1, sqrt(10) )}
                       ,I=1e5
)
y_trace4 = response4[[1]]; accpt_4 = response4[[2]]

plot(density(y_trace1))
lines(density(y_trace4), col='red')
最後，我們越來越接近我們正在尋找的東西。從理論上講，如果我們讓採樣器運行足夠長的時間，我們可以從任何這些提議分佈中得到一個有代表性的樣本，但是隨機遊走非常快地產生了一個可用的樣本，我們必須利用我們對後驗如何假設的知識尋求調整固定的採樣分佈以產生可用的結果（說實話，我們還沒有完全做到這一點y_trace4）。

稍後我將嘗試以大都會黑斯廷斯的示例進行更新。您應該能夠相當容易地看到如何修改上述代碼以生成 Metropolis Hastings 算法（提示：您只需將補充比率添加到logR計算中）。