MCMC 是Markov Chain Monte Carlo 的簡稱，但在傳統模擬中有一個很重要的假設是樣本是獨立的（independent samples），這一點在貝葉斯統計尤其是高緯度的模型中很難做到。所以MCMC的目的就是運用蒙特卡洛模擬出一個馬可鏈（Markov chain)。
如今，概率建模風靡一時，但是當我第一次了解它時，總有一件事情困擾我。 許多貝葉斯建模方法都需要計算積分，而我看到的任何工作示例似乎都使用高斯或伯努利分布，原因很簡單如果您嘗試使用比這更複雜的方法，它將成為分析的噩夢。 將貝葉斯模型限制在「表現良好」的分布的小子集中，可能會極大地阻礙你對問題建模的能力，所以我們必須找到克服這一限制的方法。
蒙特卡洛近似
如果我不想分析計算某個討厭的積分怎麼辦?可以使用蒙特卡洛近似。
我們知道，我們可以通過使用目標分布的樣本值計算期望通過使用目標分布的樣本值計算樣本均值。為什麼重要?那麼，期望是什麼呢?
連續隨機變量的期望。同樣的過程也適用於離散的情況，只要改變求和的積分。
這種估計積分的方法由中心極限定理提供了一些很好的保證。首先，這是期望的無偏估計，其次，我們可以計算估計的方差。
使用蒙特卡羅樣本計算積分是非常好的，但是我們如何從目標分布中抽取樣本呢?繪製高斯或均勻樣本很容易，但np.random會讓你失望。
畫樣本最簡單的方法是使用逆CDF方法但這依賴於獲得逆CDF函數它通常沒有一個很好的解析形式只對一維隨機變量有意義。
Metropolis算法是許多馬爾可夫鏈蒙特卡洛（MCMC）採樣方法的組成部分之一。 當您可以訪問的只是目標分布的pdf時，它使我們能夠繪製樣本。
MCMC方法需要注意的是我們不再採用獨立樣本所以我們不能保證估計的方差如何隨著樣本數量的增加而減少。如果樣本是獨立的，中心極限定理告訴我們估計值的方差將與樣本數量(N)成反比地減少。對於MCMC，我們可以通過將樣本數量從N調整為Neff來對其貼現。Neff(幾乎)總是小於N，與鏈中樣本的相關性有關。
Metropolis採樣
Metropolis算法的步驟如下：
1.從目標分布域或先前分布的域中均勻採樣起點。 2.在那時pdf。 3.根據一些狀態轉換函數，從當前位置走一步，為新樣本提出建議。 4.計算新的pdf值。 5.計算新pdf /舊pdf的值。 6.如果比率大於1，請接受該步驟。 7.如果比率小於1： 1.採樣一個統一的隨機變量。 2.如果比率大於隨機樣本，請接受該步驟。 8.否則拒絕該建議，將當前職位作為示例添加並採取新的步驟。
請注意，在5–8中描述的過程等效於以概率為min（1，p（new）/ p（old））的伯努利概率接受樣本，記住這個後面會用到……
Metropolis為何起作用？
對於任何MCMC方法，我們都希望確保一種稱為詳細平衡或可逆性的屬性。
詳細平衡
如果π滿足，則π是馬爾可夫鏈（1）的平穩分布。 我們可以通過對等式兩邊求和來證明這一點。 如果我們可以保證詳細的平衡，那麼我們也知道我們正在從馬爾可夫鏈的固定分布中取樣，我們將其作為目標分布。
詳細平衡的思想是，由於每次轉移的概率「質量」從狀態i到狀態j的轉移與從狀態j到狀態i的轉移是相同的，因此在鏈的每次轉移之後，我們保持在平穩分布 。
現在，讓我們展示一下Metropolis算法如何滿足這一條件……
為了找到滿足詳細平衡的p（i，j），我們首先提出任意的「跳躍」概率q（i，j），然後僅通過接受概率為α（的「跳躍」來獲得p（i，j）。 i，j）。 當「跳轉」被拒絕時，狀態保持j = i。 這種「接受」的想法並不是Metropolis算法獨有的，它存在於MCMC的大多數變體中。
跳躍概率推導
這取決於α是有效概率分布。 因此，α的有效形式為：
Metropolis-Hastings 跳躍概率
如果跳躍概率是對稱的，則可以簡化為：
否則，它可以保留為完整形式，稱為Metropolis-Hasting MCMC。
現在我們可以保證詳細的平衡，我們可以讓馬爾可夫鏈式接管。 如果馬爾可夫鏈是遍歷的（所有狀態都是不可約的），那麼在某個時候，該鏈將到達平穩分布，並且我們能夠從目標分布中獲取樣本。
你可能還注意到，由於alpha是π（j）/π（i）的函數。 這樣的結果意味著目標分布無需標準化。 當將Metropolis用於難以計算證據項的貝葉斯後驗估計時，這特別有用。
實現的注意事項
Metropolis算法的通用版本稱為「隨機行走Metropolis」，其中建議的狀態為當前狀態，再加上均值為零且協方差矩陣為σI的多元高斯。σ應選擇為足夠使得足夠多的樣本被拒絕大。 這是為了確保對目標分布進行良好的探索。
要注意的第二件事是老化的概念。 在馬爾可夫鏈到達平穩分布之前採集的樣本應刪除，因為它們在鏈收斂之前不能代表目標分布。 確定要刪除的樣本數量有些困難，但通常，要刪除的樣本為10–20％（或有效的樣本為10–100）。
Numpy代碼實現
def metropolis(pi, dims, n_samples, burn_in=0.1, var=1): theta_ = np.random.randn(dims)*var samples = [] while len(samples) < n_samples: theta = theta_ + np.random.randn(dims)*varratio = pi(theta)/pi(theta_) if np.random.rand(1) < ratio: sample = theta theta_ = theta samples.append(sample) else: sample = theta_ samples.append(sample) samples = np.array(samples) return samples[int(samples*burn_in):,:]
我們可以看到這在兩個高斯函數的和上的表現(注意這是一個非正態分布)。
from scipy.stats import multivariate_normal def make_pdf(mean1, mean2, cov1, cov2): pdf1 = multivariate_normal(mean1, cov1) pdf2 = multivariate_normal(mean2, cov2) def pdf(x): return pdf1.pdf(x) * pdf2.pdf(x) return pdfpdf = make_pdf( [3, 3], [-1, -1], np.array([[1,0.1],[0.1,1]], dtype=float), np.array([[1,0.1],[0.1,1]], dtype=float), ) samples = metropolis(pdf, 2, 1_000, 0.1)
上面的gif顯示了算法是如何遍歷分布的，偶爾會在分布的兩種不同模式之間跳轉。注意，這也突出了metropolis算法的一個弱點，它處理相對較差的多模型分布。這是由於要探索一種新模式，步驟必須足夠大，以便從一種模式希望到另一種模式。這要麼需要一個大的步長，要麼需要一個模態緊密分布。修改如哈密頓量MCMC可以幫助解決這一問題，但一般來說，這是大多數MCMC方法的一個問題。
論文地址：arxiv/1504.01896.pdf
作者:Alexander Bailey
deephub翻譯組