一、引入
拒絕采樣，重要性采樣的效率在高維空間很低，隨維度增長其難度也指數型增長，主要適用於一維的采樣。對於二維以上可以用馬氏鏈。馬爾可夫鏈蒙特卡洛采樣方法就是在高維空間采樣的方法。
馬爾可夫鏈就是滿足馬爾可夫假設的一組狀態序列$\left \{ x_{t} \right \}= ...,x_{t-2}, x_{t-1}, x_{t}, x_{t+1},...$，其中假設某一時刻狀態$x_{t}$發生狀態轉移的概率只依賴於它的前一個狀態$x_{t-1}$，這個就是馬爾可夫假設：
$P(X_{t+1}\mid X_{t-2}, X_{t-1}, X_{t}, X_{t+1} )=P(X_{t+1}\mid X_{t})$
只要我們能知道系統中任意兩個狀態之間的轉移概率，整個馬爾可夫模型就知道了，定義轉移概率為：
$P_{ij}=P(X_{t+1}=j\mid X_{t}=i)$
轉移概率衡量的是兩個狀態之間的轉移幾率，與時刻無關，僅涉及相鄰的兩個狀態，如果系統中的狀態有T種，那么轉移概率可以構成一個T×T的轉移矩陣P，馬爾可夫鏈有收斂性質，就是說從任意一個初始分布出發，經過多次轉移后，得到的分布是平穩的，不再變化的分布q，這個平穩分布q與初始分布無關，只與狀態轉移矩陣P有關。
接下來用一段代碼說明這個問題：
matrix1=np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]],dtype=float) vector1=np.matrix([[0.3,0.4,0.3]],dtype=float) for i in range(100): vector1=vector1*matrix1 print('current round:',i+1) print(vector1)
輸出最開始定義的matrix1和vector1，分別是轉移概率和初始概率分布：
初始概率分布vector1表示的是在t=0時刻，$P(X_{0}=0)=0.3$，$P(X_{0}=1)=0.4$，$P(X_{0}=2)=0.3$，描述的是起初時刻，狀態的分布情況，這里默認有三種狀態，分別是1,2,3。
轉移矩陣matrix1的位置11上的元素0.9表示的是：當前一時刻的狀態$X_{t-1}$取1時，下一時刻的狀態$X_{t}$取1的概率為0.9；
位置12上的元素0.075表示的是：當前一時刻的狀態$X_{t-1}$取1時，下一時刻的狀態$X_{t}$取2的概率為0.075；
可以用vector1*matrix1用矩陣乘法公式來計算一下，會得到一個形式和vector1一樣的矩陣，有三個元素，每個元素表示的就是該時刻的狀態分布情況：
根據矩陣乘法公式來計算一下下一時刻的狀態分布的第一個元素為：（0.3*0.9+0.4*0.15+0.3*0.25），可以看出，它表示的是【不管上一時刻是取了哪個值，下一時刻狀態為1的概率】，這個元素的計算考慮了上一個時刻的狀態分布概率，比如0.3*0.9，這個乘積的含義就是上一時刻狀態為1的概率是0.3，在上一時刻狀態為1的條件下下一時刻狀態仍為1的概率是0.9，注意這里的0.9是個條件概率；相乘就可以得到上一時刻狀態為1，且下一時刻狀態為1的概率，這里注意這是個聯合概率分布。
輸出一下整體的代碼結果：
我們發現最后得到的分布基本是平穩不變的，即：狀態取0和3的概率是0.625，狀態取1的概率是0.3125。我們可以改變初始vector1，但是最終收斂的平穩分布依舊是不變的。這個就是馬氏鏈的收斂。
我們的目標是希望從平穩分布q中進行采樣。現在的問題如下：
- 什么樣的馬氏鏈可以收斂到一個平穩分布q？所有的馬氏鏈都收斂嗎？這個在截圖中已有。
- 既然q只和轉移矩陣P有關，我們要如何確定馬氏鏈中的轉移矩陣P，使得平穩分布是我們想要的q？在三中說明
- 我們具體要怎么從平穩分布中采樣？這個在二中大體說明，在三中具體說明
二、怎么基於馬爾可夫鏈采樣？
假設我們有初始概率分布為$\pi _{0}\left ( x \right )$，基於$\pi _{0}\left ( x \right )$采樣得到$x_{0}$；
然后基於條件概率分布$P\left ( x\mid x_{0} \right )$采樣得到$x_{1}$；
基於條件概率分布$P\left ( x\mid x_{1} \right )$采樣得到$x_{2}$；
當這個過程進行到第n個時刻時，假設此時達到平穩分布，即條件概率不再變化，如果我們的目標是要采m個樣本的話，那么就從第n個時刻開始，再根據平穩分布$q$采m次得到m個樣本，再做蒙特卡洛求和即可得到我們的目標期望。
這節只是為了引出下一節。在后面還會提到更具體的算法來說明如何從平穩分布里采樣。
但是往往，我們只知道我們要從一個分布$q$中采樣，這個$q$很難采我們才去找馬爾可夫鏈，因為馬爾可夫鏈進行到最后所得到的樣本服從我們的目標分布$q$。但是$q$是多少主要取決於$P$，我們要怎么去確定$P$？這就是下一節會解決的第二個問題。
這里還要說明的就是，我們知道目標分布$q$的表達式，但這並不意味這我們可以直接從中采樣，因為很難采，所以引入馬爾可夫鏈，將直接采樣方式轉換為較為簡單的條件概率分布采樣。
三、馬爾可夫鏈的細致平穩條件
上面代碼部分說明了馬氏鏈不斷進行更新，最后會得到一個平穩分布，但是那是基於大量的計算才得到的結果。到底怎么樣才算是一個平穩分布，現在給出定義：
如果非周期馬爾可夫鏈的狀態轉移矩陣P和概率分布$\pi(x)$滿足$\pi(i)P_{ij}=\pi(j)P_{ji}$，那么就稱\pi(x)是狀態轉移矩陣P的平穩分布，再進一步對該式兩邊對狀態i進行求和可得：
$\sum_{i}^{\infty }\pi(i)P_{ij}=\sum_{i}^{\infty }\pi(j)P_{ji}=\pi(j)\sum_{i}^{\infty }P_{ji}=\pi(j)$
上式證明了滿足$\pi(i)P_{ij}=\pi(j)P_{ji}$的話，的確可以收斂。這個條件稱為馬爾可夫鏈的細致平穩條件，滿足了這個條件的分布才能稱為平穩分布。這個式子較為簡潔的搭建了平穩分布和狀態轉移矩陣之間的關系。可以用來解決我們上面提到的第二個問題。
首先規定一個目標分布$\pi(x)$是我們希望得到的平穩分布，也就是很難直接采樣的分布$q$，和一個馬爾可夫鏈狀態轉移矩陣$Q$，顯然最開始它們是不滿足細致平穩條件的：
$\pi(i)Q_{ij}\neq \pi(j)Q_{ji}$
我們引入一個$\alpha (ij)$，使得上式成立：
$\pi(i)Q_{ij}\alpha (ij)= \pi(j)Q_{ji}\alpha (ji)$
其中有：
$\alpha (ij)= \pi(j)Q_{ji}$
$\alpha (ji)=\pi(i)Q_{ij}$
所以我們得出一個結論：目標平穩分布對應的概率轉移矩陣$P$，可以由任意一個狀態轉移矩陣$Q$，乘上一個接受率$\alpha (ij)$得到；也就是說，想要得到的$P$，可以通過不正確的$Q$以一定的概率獲得。所以我們在二、的基礎上進一步細化我們馬爾可夫鏈的采樣過程如下：
1）初始擁有：任意轉移矩陣Q，平穩分布$\pi(x)$，轉移次數$n_{1}$，所需樣本數$n_{2}$
2）任意采樣初始狀態值$x_{0}$
3）$for t=0 to n_{1}+n_{2}-1$：
a）從條件概率分布$Q\left ( x\mid x_{t} \right )$中采樣得到樣本$x_{*}$
b）采樣$u\sim U[0,1]$
c）如果u< \alpha (x_{t},x_{*})=\pi(x_{*})Q(x_{*},x_{t})，就接受轉移，即采樣$x_{*}$，否則不接受，有$x_{t+1}=x_{t}$
為什么這里要采樣一個均勻分布u樣本點，來和接受率進行比較呢？因為接受率反應了一件事——原來的轉移矩陣Q不符合我們的需求，它太大了，有一些樣本不符合我們的所需，需要丟掉，但是怎么丟，就是看接受率，只有很少一部分的樣本可以被我們選到。采樣均勻分布體現了一個任意性，每個樣本的接受率都是均勻的，意思是不存在有100個樣本的接受率是0.9，而只有1個樣本的接受率是0.1，均勻分布就是避免了這種情況，說明我們的樣本集非常完整。就像如果要統計某病在人群中的發病率，我們最好控制變量，最好我們的樣本集男女數目一樣多，全國各省人的比例也和地區人數呈正比。
還有這里的轉移次數$n_{1}$要怎么確定的問題，轉移次數其實可以設置得大一些，因為我們采的樣本是從第$n_{1}+1$次開始采的，前面都在使馬爾科夫鏈達到平穩，只有平穩了我們才能開始收集樣本。
這個方法也有缺點，就是如果我們的接受率很小$\alpha (x_{t},x_{*})$，使得我們的大部分采樣點都被拒絕了，這樣的話，采樣效率很低，也就是說會出現很多$x_{t+1}=x_{t}$的情況，采樣樣本中有很多相同的點，這不利於我們拿這些樣本來做蒙特卡洛求和，因為很有可能這些樣本對應的函數值很小，對整體起不到很大的作用。