文章目錄

• 什么是采樣
• 為什么要采樣
• 采樣的方法
• • 蒙特卡羅采樣
  • 重要采樣
  • 馬爾科夫鏈蒙特卡羅方法（MCMC）
  • • 馬氏過程相關概念和定理
    • MCMC算法
    • Metropolis-Hastings算法
    • MCMC方法的問題
  • 吉布斯采樣（Gibbs）
• 參考文獻

什么是采樣

我們知道了一個變量的分布，要生成一批樣本服從這個分布，這個過程就叫采樣。比如在信號處理領域，采樣是將信號從連續時間域上的模擬信號轉換到離散時間域上的離散信號的過程。

為什么要采樣

有許多原因需要我們從某個分布中采樣。比如獲取離散信號，比如數據挖掘中進行模型匹配，比如機器學習中估計某個函數在某個分布上的期望，比如以較小的代價近似許多項的和或者某個積分。
在許多情況下，采樣的目標就是抽樣，比如對大數據集中獲取合適的子集。

采樣的方法

蒙特卡羅采樣

蒙特卡羅方法主要是找到滿足分布的隨機變量，對該隨機變量隨意取值生成若干樣本。當樣本足夠多時，樣本分布近似于目標分布。
當我們無法精確計算和或積分時，可以采用蒙特卡羅方法來近似。這種想法是把和或積分視作某個分布下的期望，然后通過估計對應的平均值來近似這個期望。令
$$s=\sum _{x}p(x)f(x)=E_p[f(x)]$$
或者
$$s=\int p(x)f(x)dx=E_p[f(x)]$$
其中$p(x)$為概率分布函數(求和)或者概率密度函數(積分)
我們可以通過從$p$中抽取n個樣本$x^{(1)},...,x^{(n)}$來近似s并得到一個經驗平均值
$$s_n=\frac{1}{n}\sum _{i=1}^{n}f(x^{(i)})$$
顯然這是一個無偏估計，且由大數定律，$\underset{n\to \infty }{\lim }{s}_n=s$.
由中心極限定理有$s_n～\mathcal{N}(s,\frac{Var[f(x)]}{n})$，由此可以估計$s_n$的置信區間。

重要采樣

蒙特卡羅采樣中假設了可以從基準分布$p(x)$中輕易采樣，當無法輕易采樣時，可以換一個能夠容易采樣的分布，如上例中，$p(x)f(x)=q(x)\frac{p(x)f(x)}{q(x)}$，由此可以從$q$分布中采樣，然后估計$\frac{pf}{q}$的期望。這樣的轉化顯然不唯一，因此我們可以尋找最優的分布$q^{?}$。這樣的分布可以容易的推導出來。令
$$s_q=\frac{1}{n}\sum _{i=1,x^{(i)}～q}^n\frac{p(x^{(i)})f(x^{(i)})}{q(x^{(i)})}$$
計算得：
$$\mathbb{E}[s_q]=s,Var[s_q]=Var[\frac{p(x)f(x)}{q(x)}]/n$$
可以發現估計的期望與$q$分布無關。我們希望方差取最小值，則$q$需要滿足：
$$q^{?}(x)=\frac{p(x)|f(x)|}{Z}$$
這里$Z$是歸一化常數，此時$Var[s_q]=0$.

一個好的$q$分布的選擇可以顯著地提高蒙特卡羅估計的效率，而一個糟糕的$q$分布則會使效率更糟糕。q分布通常會取一些簡單常用的分布使得我們能夠從$q$分布中容易的采樣，但當$x$是高維數據時，$q$分布的簡單性使得它很難和$p$或$p|f|$相匹配。

馬爾科夫鏈蒙特卡羅方法（MCMC）

馬爾科夫蒙特卡羅方法的基本思想是構造一個馬氏鏈$\{X_n:n?0\}$，使目標分布作為該馬氏鏈的不變分布，再利用馬氏鏈的極限理論（遍歷定理），用$X_n$代替目標分布。

馬氏過程相關概念和定理

定義1 轉移矩陣
設狀態空間$E$上的矩陣$P=(p(i,j):i,j\in E)$為轉移矩陣，滿足一下性質：
(1) 對任何$i,j\in E$，有$p(i,j)?0$
(2) 對任何$i\in E$，有$\sum _{j\in E}p(i,j)=1$

定義2 馬氏鏈
對于隨機過程$X=\{X_0,X_1,X_2,...,X_n,...\}$。我們稱$X$是以$P$為轉移矩陣的馬氏鏈，當且僅當$P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)$
由此可以得到$P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)=P(X_{n+1}=j|X_n=i_n)$
特別地，利用轉移矩陣和初始分布可以給出馬氏鏈的有限維分布的表示。

定義3 平穩分布（不變分布）
對于以$P$為轉移矩陣的馬氏鏈$X$，如果概率分布$\eta =({\eta }_i:i\in E)$滿足$\sum _{i\in E}{\eta }_i{p}_{i,j}={\eta }_j,j\in E$，則稱$\eta$為$P$或$X$的平穩分布（不變分布）。
定義4 可逆分布
對于以$P$為轉移矩陣的馬氏鏈$X$，如果概率分布$\eta =({\eta }_i:i\in E)$滿足${\eta }_i{p}_{i,j}={\eta }_j{p}_{j,i},i,j\in E$，則稱$\eta$為$P$或$X$的可逆分布。
顯然可逆分布一定是平穩分布。

如果馬氏鏈X以$P$為轉移矩陣，以平穩分布$\eta$為初始分布，則稱馬氏鏈$X$是平穩的。
如果馬氏鏈X以$P$為轉移矩陣，以可逆分布$\eta$為初始分布，則稱馬氏鏈$X$是可逆的，也叫細致平衡的。

定理1（弱遍歷定理）
設轉移矩陣為$P$，令$L^{(n)}=\frac{1}{n}\sum _{n=0}^n{P}^{(n)}$。則$\underset{n\to \infty }{\lim }{L}^{(n)}=L$存在且滿足$PL=LP=L=L^2$

定理2
設$P$是馬氏鏈的轉移矩陣，假定$P$不可約正常返，則$P$有唯一不變分布，并且轉移概率的平均極限分布是馬爾可夫鏈的平穩分布。

事實上
$$\underset{n\to \infty }{\lim }{L}^{(n)}=L=?\begin{array}{c}\pi \\ \pi \\ .\\ .\\ .\end{array}?$$
$\pi$為行向量且任意元素都大于0，此時$\pi$為$P$的唯一不變分布。

定理3
若$P$不可約，正常返，非周期，則有$\underset{n\to \infty }{\lim }=P^{(n)}=L$

MCMC算法

有了這些理論就可以給出馬爾科夫蒙特卡羅法的采樣過程了。先假設存在轉移矩陣$P$，目標分布就是$P$的不變分布。
先從初始分布從采一個樣本$x^{(0)}$，這個初始分布可以是狀態空間上的任意分布。接著從利用轉移矩陣得到在$x^{(0)}$條件下的分布$\mathbf{P}(?|X_0=x^{(0)})$，在此分布下采樣得到$x^{(1)}$，不斷重復上述過程，可以發現每個時刻在這個馬爾可夫鏈上進行隨機游走一次，就可以得到一個樣本。

# 算法1：生成一條馬氏鏈軌跡 Proceduce MCMC_Sample(E #狀態空間P_0 #初始分布P #馬爾科夫轉移模型n #執行步數 ) 1 Sample x^(0)~P_0 2 for t=1,2,...,n 3 Sample x^(t)~P(·|x^(t-1)) 4 return x^(0),x^(1),...,x^(n)

在算法執行$m(<n)$步后，可以認為得到的樣本集合$\{x^{(m+1)},x^{(m+2)},...x^{(n)}\}$為所需要的樣本。前$m$步稱為混合期。
現在我們可以得到對$s$的估計$s_M=\frac{1}{n?m}\sum _{i=m+1}^{n}f(x^{(i)})$

Metropolis-Hastings算法

在MCMC算法中，我們假設了已經有這樣的轉移矩陣$P$，但事實上這樣的矩陣很難直接找到。而Metropolis-Hastings算法就是用一種簡單的方法構造出來轉移矩陣$P$，從而能夠應用MCMC算法采樣。

注意到細致平衡的分布也是不變分布，由此可以構造一個符合條件的可逆過程的轉移矩陣。設目標分布為$\pi$，將目標分布設為轉移矩陣$P$的不變分布。首先令$Q$是任意一個取正整數的特定不可約馬爾科夫轉移概率矩陣，以$q(i,j)$表示$Q$的$i$行$j$列的元素。設目標分布為$\pi$現按照細致平衡的條件，對任意$i=j$，有
$$\pi (i)q(i,j)a(i,j)=\pi (j)q(j,i)a(j,i)$$
其中$a(i,j)$為配比出來的系數，于是可令$a(i,j)=\min (\frac{\pi (i)q(i,j)}{\pi (j)q(j,i)},1)$
所以可設$p(i,j)=q(i,j)a(i,j)$，為保持轉移矩陣的性質，得到$p(i,i)=q(i,i)+\sum _{k=i}q(i,k)(1?a(i,k))$，最后設$P=(p(i,j))$，則$P$即為所需的轉移矩陣。

# 算法2：M-H算法 Proceduce MCMC_Sample(E, P_0, Q, a, n) 1 Sample x^(0)~P_0 2 for t=1,2,...,n 3 Sample x*~Q(·|x^(t-1)) 4 Sample u~U[0,1] 5 if u<a(x^(t-1),x*) then 6 x^(t) = x* 7 else 8 x(t) = x(t-1) 9 return x^(0),x^(1),...,x^(n)

MCMC方法的問題

采集的相鄰樣本可能含有相關性
一種解決辦法是每隔n個樣本返回一個樣本，但這樣不僅提高了采樣時的計算代價，而且也沒有完全解除掉相關性；另一種解決辦法是同時并行多個馬氏鏈，每隔馬氏鏈只采集一個樣本。這兩種方法是兩個極端，更多從業者選擇馬氏鏈和小批量中的樣本數接近，然后再用這些馬氏鏈抽取所需要的樣本。

無法預測混合期
檢測一個馬爾科夫鏈是否達到了平衡很困難，只能粗略的估計這段時間是否是足夠的。由此可見，此算法的關鍵就是能否針對具體的$\pi$找到一個能夠快速收斂的馬氏鏈，收斂速度的快慢決定了算法的優劣。
由不變分布方程$\pi P=\pi$以及不變分布的唯一性知，$P$的特征值中有且只有一個$1$，且$\pi$是對應的唯一特征向量。注意到$P$是不可約的，則由Perron-Frobenius定理可以保證這個矩陣最大特征值為1.這就導致了其他特征值都隨時間變化趨向于0，那么第二大特征值將決定了收斂速度，也就是混合期的大小。
注： 在M-H算法中，預設分布的標準差也影響著混合時間。如果標準差過小，則采集的樣本多集中在概率大的周圍；若標準差過大，則被拒絕的概率將變大，混合時間將變長。

吉布斯采樣（Gibbs）

在計算過程中常常會遇到復雜的聯合分布，這增加了計算的復雜度。事實上不同分量之間可能存在相關性，利用概率圖模型，將聯合分布分解成多個條件分布，可以減少狀態空間的參數。在采樣過程中，將高維總體的采樣，化為一系列一維分布的取樣，這就是Gibbs采樣法的要義。

# 算法3：Gibbs算法 Proceduce MCMC_Sample(E, P_0, P, T) 1 Sample x^(0)~P_0 2 for t=1,2,...,T 3 Sample x_1^(t) ~ P(x_1|x_2^(t-1), x_3^(t-1),...,x_n^(t-1)) 4 Sample x_2^(t) ~ P(x_2|x_1^(t), x_3^(t-1),...,x_n^(t-1)) 5 ... 6 Sample x_j^(t) ~ P(x_j|x_1^(t),..., x_(j-1)^(t),x_(j+1)^(t-1),...,x_n^(t-1)) 7 ... 8 Sample x_n^(t) ~ P(x_n|x_1^(t), x_2^(t),...,x_(n-1)^(t)) 9 return x^(0),x^(1),...,x^(n)

參考文獻

古德費洛, 本吉奧, 庫維爾, 趙申劍, 黎彧君, and 符天凡. 深度學習. 北京: 人民郵電出版社, 2017.

錢敏平. 應用隨機過程. 北京: 高等教育出版社, 2011.

羅斯, and 龔光魯. 應用隨機過程 概率模型導論. 北京: 人民郵電出版社, 2011.

龔光魯. 應用隨機過程教程 及在算法和智能計算中的隨機模型. 北京: 清華大學出版社, 2004.

Andrieu, Christophe, De Freitas, Nando, Doucet, Arnaud, and Jordan, Michael I. “An Introduction to MCMC for Machine Learning.” Machine Learning 50.1 (2003): 5-43.

科勒, 弗里德曼, 王飛躍, and 韓素青. 概率圖模型 原理與技術. 北京: 清華大學出版社, 2015.

Perron, Oskar. “Zur Theorie Der Matrices.” Mathematische Annalen 64.2 (1907): 248-63.

總結

以上是生活随笔為你收集整理的马尔科夫蒙特卡罗方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。