采樣之MCMC采樣和M-H采樣

在采樣之馬爾科夫鏈中我們講到給定一個概率平穩分布π, 很難直接找到對應的馬爾科夫鏈狀態轉移矩陣P。而只要解決這個問題，我們就可以找到一種通用的概率分布采樣方法，進而用於蒙特卡羅模擬。本篇我們就討論解決這個問題的辦法：MCMC采樣和它的易用版M-H采樣

1.馬爾科夫鏈的細致平穩條件

2. MCMC采樣

假設我們已經有一個轉移矩陣Q(對應元素為q(i,j)), 把以上的過程整理一下，我們就得到了如下的用於采樣概率分布p(x)的算法。

3. M-H采樣

4. M-H采樣實例

為了更容易理解，這里給出一個M-H采樣的實例。

　　　　在例子里，我們的目標平穩分布是一個均值3，標准差2的正態分布，而選擇的馬爾可夫鏈狀態轉移矩陣Q(i,j)的條件轉移概率是以i為均值,方差1的正態分布在位置j的值。這個例子僅僅用來讓大家加深對M-H采樣過程的理解。畢竟一個普通的一維正態分布用不着去用M-H采樣來獲得樣本。

　　　　代碼如下：

# coding: utf-8
import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt

# scipy.stats.norm.pdf(x, loc=0, scale=1)
# 輸入x，返回概率密度函數;loc代表了均值,scale代表標准差
def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    # rvs產生服從正太分布的一個樣本，對隨機變量進行隨機取值
    # rvs可以通過size參數指定輸出的數組大小，即如果size=1 則產生一個樣本值，如果size=2，則產生兩個樣本值。
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]

# scatter是制作x與y的散點圖
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7)
plt.show()

輸出的圖中可以看到采樣值的分布與真實的分布之間的關系如下:

 

 參考：https://www.cnblogs.com/xbinworld/p/4266146.html

           http://www.cnblogs.com/pinard/p/6638955.html