Markov chain
基本概念
當随機變量序列\(\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}\}\)的聯合分布滿足：
\(p(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n})=p(\mathbf{x}_{1})\prod_{t=2}^{n}p(\mathbf{x}_{t}|\mathbf{x}_{t-1})\)
時，稱其為Markov model或者Markov chain。特點是\(\mathbf{x}_{t}\)包含了所有過去的資訊，下一随機變量
\(\mathbf{x}_{t+1}\)的分布完全由\(\mathbf{x}_{t}\)決定，\(p(\mathbf{x}_{t}|\mathbf{x}_{1:t-1})=p(\mathbf{x}_{t}|\mathbf{x}_{t-1})\)。
\(\mathbf{x}_{t}\)的邊緣分布為：
\(\mathbf{x}_{t}=\int_{\mathbf{x}_{1:t-1}}p(\mathbf{x}_{1})\prod_{t=2}^{n}p(\mathbf{x}_{t}|\mathbf{x}_{t-1})d\mathbf{x}_{1}d\mathbf{x}_{2}...d\mathbf{x}_{t-1} \)
當\(\mathbf{x}_{t}\)的值域為有限集合\(\{1,2,\cdots,K\}\)時，這種離散的馬爾可夫鍊可以用\(K\times K\)躍遷矩陣\(A\)表示，其中\(A_{ij}=p(\mathbf{x}_{t}=j|\mathbf{x}_{t-1}=i)\)表示從狀态\(i\)躍遷到狀态\(j\)的機率。
另外，用\(A(n)\)表示\(n\)步躍遷，即\(A_{ij(n)}=p(\mathbf{x}_{t+n}=j|\mathbf{x}_{t}=i)\)。
對于離散情況，利用躍遷矩陣\(A\)，可以将\(\mathbf{x}_{t}\)的邊緣分布表示：
\(\mathbf{x}_{t}=\int_{\mathbf{x}_{1:t-1}}p(\mathbf{x}_{1})\prod_{t=2}^{n}p(\mathbf{x}_{t}|\mathbf{x}_{t-1})d\mathbf{x}_{1}d\mathbf{x}_{2}...d\mathbf{x}_{t-1}=p(\mathbf{X}_{1})A\times A...\times A=p(\mathbf{X}_{1})A^{n-1}\)
上式利用了關系\(p(\mathbf{x}_{t})=\int p(\mathbf{x}_{t}|\mathbf{x}_{t-1})d\mathbf{x}_{t-1}=p(\mathbf{x}_{t-1})A\)
Markov chain的性質
MCMC主要利用的是滿足某些條件的Markov chain具有stationary distribution的性質。這一部分可以參考文獻 [4] 的17.2.3節(簡單易懂)以及文獻[3]的第6章(更詳細)。
MCMC方法用到的一個定理(文獻[4], Theorem 17.2.3)：
If a Markov chain with transition matrix \(A\) is /regular/ and satisfies /detailed balance/ wrt distribution \(\pi\), then \(\pi\) is a stationary distribution of the chain.
MCMC
MCMC相比 Reject Sampling and Importance Sampling的主要優點是MCMC在高維的情況下的效果更好。MCMC的基本思路是，已知待采樣的分布\(\pi(\mathbf{x})\)，通過構造躍遷核\(A(\mathbf{x},\mathbf{x}')\)，得到一個馬爾可夫鍊\(M\)。\(M\)具有規則性(regular)并且對\(\pi(\mathbf{x})\)滿足detailed balance條件。這樣，\(M\)具有穩态分布\(\pi(\mathbf{x})\)。對\(M\)依次進行采樣，得到采樣點序列\(\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n},\cdots\}\)。當采樣點數\(n\)足夠大時，可以認為\(p(\mathbf{x}_{m>n})\)已經進入了穩态分布，進而\(\mathbf{x}_{m>n}\)就是從穩态分布\(\pi(\mathbf{x})\)得到的采樣。注意這樣得到的采樣點是關聯的(通常非常相關)，這一點和reject sampling不同。
上述“對\(M\)進行采樣”就是對分布\(A(\mathbf{x}=\mathbf{x}_{0},\mathbf{x}')\)進行采樣，其中\(\mathbf{x}_{0}\)是目前已經得到的采樣點。
在離散情況下，\(A(\mathbf{x},\mathbf{x}')\)退化成一個矩陣，分布\(A(\mathbf{x}=\mathbf{x}_{0}=j,x')\)就是矩陣\(A\)的第\(j\)行定義的機率分布\(p(\mathbf{x}_{t+1}|\mathbf{x}_{t}=j)\)
在具體實作的時候，關于\(A\)的構造和采樣可能很難實作。我們可以采用類似于Reject Sampling的方法：選擇一個很簡單的分布(proposal distribution)\(T=q(\mathbf{x}_{t+1}|\mathbf{x}_{t})\)進行采樣，然後再對采樣點進行判斷，以機率\(r\)接受該采樣。這樣一來，\(A=Tr\)由兩部分組成，可以設計\(r\)的形式(類似Reject Sampling，\(r\)一般由\(\pi\)和\(T\)同時決定)，使\(A\)滿足上文所述的兩個條件。典型的此類算法有Metropolis算法和Metropolis-Hastings算法。
Metropolis algorithm的步驟如下：
1選擇proposal distribution \(T=q(\mathbf{x}_{t+1}|\mathbf{x}_{t})\)。\(T\)的常見選擇可以參考random-walk Metropolis、Metropolized independent sampler等等
令\(t=0\)，任選一個初始采樣點\(\mathbf{x}_{0}\)，比如從均勻分布中采樣得到一個樣本
對分布\(T\)進行采樣，得到候選采樣點\(\bar{\mathbf{x}}\)
計算\(r=min\{1,\frac{p(\bar{\mathbf{x}})}{p(\mathbf{x}_{t})}\}\)
以機率\(r\)接受候選采樣點\(\bar{\mathbf{x}}\)。如果接受，則\(\mathbf{x}_{t+1}=\bar{\mathbf{x}}\);如果拒絕，則\(\mathbf{x}_{t+1}=\bar{\mathbf{x}}\)
重複3-5直到得到所需采樣點數
Metropolis算法要求proposal distribution滿足對稱性，進而增加了proposal distribution的選取難度。Hastings對Metropolis算法做了改進，取消了proposal distribution必須對稱的要求，該算法即Metropolis-Hastings，具體細節可以參見[3]或[1]。
Gibbs算法
Gibbs算法可以認為是一類特殊的MCMC算法，在統計實體中被廣泛應用。從形式上來看，Gibbs算法把對一個\(N\)維變量的采樣分解為\(N\)個一維的采樣，即每次對變量的一個坐标進行采樣，得到一個采樣點需要進行\(N\)次采樣過程。可以證明Gibbs算法每次采樣的接受率都是1，即從來不會拒絕采樣點。
Gibbs算法的突出優點是它的每一步都是在進行一維采樣，這一般是非常簡單的。比如在機率圖模型中的使用partical based inference方法計算某些機率分布時，就可以利用Gibb算法友善的獲得采樣點。因為進行采樣時除了待采樣變量，其它變量的值都是固定的，可以很友善的求出待采樣變量的一維分布，然後進行采樣。
下面是Metropolis算法和Gibbs采樣算法的一個簡單例子，來自 Gibbs sampler
1 #-*- coding: utf-8 -*-
2
3 """
4 利用MCMC算法對兩維分布采樣5 問題描述：6 給出兩個硬币A、B，進行抛硬币實驗，其中A硬币共抛擲N_1次，正面向上次數z_1次，B硬币共抛擲N_2次，正面向上次數z_2次。7 現在估計兩個硬币正面向上的機率θ_1和θ_2，即p(θ|X)=p(θ_1, θ_2|N_1, z_1, N_2, z_2)。8 已知：9 先驗分布p(θ)是Beta分布，10 似然p(X|θ)是Bernoulli分布11 後驗分布p(θ|X)是Beta分布12 兩個硬币獨立進行實驗，是以機率可以表示成兩個硬币各自分布乘積的形式13
14 source: https://people.duke.edu/~ccc14/sta-663/MCMC.html#gibbs-sampler15
16 """
17 importsys18 importglob19 importnumpy as np20 importpandas as pd21 importscipy.stats as stats22 importmatplotlib.pyplot as plt23 from mpl_toolkits.mplot3d importAxes3D24 from functools importpartial25
26
27 defbern(theta, z, N):28 """Bernoulli likelihood with N trials and z successes."""
29 return np.clip(theta**z * (1-theta)**(N-z), 0, 1)30
31
32 defbern2(theta1, theta2, z1, z2, N1, N2):33 """Bernoulli likelihood with N trials and z successes."""
34 return bern(theta1, z1, N1) *bern(theta2, z2, N2)35
36
37 defmake_thetas(xmin, xmax, n):38 xs =np.linspace(xmin, xmax, n)39 widths = (xs[1:] - xs[:-1])/2.0
40 thetas = xs[:-1] +widths41 returnthetas42
43
44 def make_plots(X, Y, prior, likelihood, posterior, projection=None):45 fig, ax = plt.subplots(1, 3, subplot_kw=dict(46 projection=projection, aspect='equal'), figsize=(12, 3))47 if projection == '3d':48 ax[0].plot_surface(X, Y, prior, alpha=0.3, cmap=plt.cm.jet)49 ax[1].plot_surface(X, Y, likelihood, alpha=0.3, cmap=plt.cm.jet)50 ax[2].plot_surface(X, Y, posterior, alpha=0.3, cmap=plt.cm.jet)51 else:52 ax[0].contour(X, Y, prior)53 ax[1].contour(X, Y, likelihood)54 ax[2].contour(X, Y, posterior)55 ax[0].set_title('Prior')56 ax[1].set_title('Likelihood')57 ax[2].set_title('Posteior')58 plt.tight_layout()59
60
61 thetas1 = make_thetas(0, 1, 101)62 thetas2 = make_thetas(0, 1, 101)63 X, Y =np.meshgrid(thetas1, thetas2)64
65 a = 2
66 b = 3
67
68 z1 = 11
69 N1 = 14
70 z2 = 7
71 N2 = 14
72
73
74 #Analytical
75
76 prior = stats.beta(a, b).pdf(X) *stats.beta(a, b).pdf(Y)77 likelihood =bern2(X, Y, z1, z2, N1, N2)78 posterior = stats.beta(a + z1, b + N1 - z1).pdf(X) *\79 stats.beta(a + z2, b + N2 -z2).pdf(Y)80 make_plots(X, Y, prior, likelihood, posterior)81 make_plots(X, Y, prior, likelihood, posterior, projection='3d')82
83
84 #采樣算法
85 defprior(theta1, theta2):86 return stats.beta(a, b).pdf(theta1) *stats.beta(a, b).pdf(theta2)87
88
89 lik = partial(bern2, z1=z1, z2=z2, N1=N1, N2=N2)90
91
92 deftarget(theta1, theta2):93 return prior(theta1, theta2) *lik(theta1, theta2)94
95
96 theta = np.array([0.5, 0.5])97 niters = 10000
98 burnin = 500
99 sigma = np.diag([0.2, 0.2])100
101 #Metropolis采樣計算p(θ|X)
102
103 thetas = np.zeros((niters-burnin, 2), np.float)104 for i inrange(niters):105 new_theta =stats.multivariate_normal(theta, sigma).rvs()106 p = min(target(*new_theta)/target(*theta), 1)107 if np.random.rand() <108 theta="new_theta109" if i>=burnin:108>
# NOTE：這裡不應該把reject的資料放到采樣結果裡110 thetas[i-burnin] =theta111
112
113 kde =stats.gaussian_kde(thetas.T)114 XY =np.vstack([X.ravel(), Y.ravel()])115 posterior_metroplis =kde(XY).reshape(X.shape)116 make_plots(X, Y, prior(X, Y), lik(X, Y), posterior_metroplis)117 make_plots(X, Y, prior(X, Y), lik(X, Y), posterior_metroplis, projection='3d')118
119
120 #Gibbs采樣計算p(θ|X)
121
122 thetas = np.zeros((niters-burnin, 2), np.float)123 for i inrange(niters):124 #θ是二維的，是以要進行兩次采樣，由于兩個參數互相獨立，是以采樣實際上是各自進行的
125 #接受率為1
126 theta = [stats.beta(a + z1, b + N1 - z1).rvs(), theta[1]]127 theta = [theta[0], stats.beta(a + z2, b + N2 -z2).rvs()]128
129 if i >=burnin:130 thetas[i-burnin] = theta
MCMC和Importance Sampling思路的差別
這兩者都是為了解決針對複雜分布的積分問題。由于目标分布太過于複雜，不但解析形式的積分沒法計算，連對其采樣都很困難。MCMC的思路是通過構造馬爾可夫鍊的方式，得到 *目标分布* 的 *關聯* 的樣本，而Importance Sampling的思路則是轉而對一個簡單的分布( *測試分布* )進行采樣，得到 *獨立* 的樣本，然後在計算積分的時候對采樣點配置設定适當的權重，進而得到對原始積分的近似解。
MCMC的收斂性讨論
使用MCMC時要注意的重要一點就是要等到MC收斂到穩态分布之後再進行真正的采樣，之前的采樣點都棄之不用，這個過程叫做 burn in。但是在實際使用中，MCMC到底何時進入穩态是難以判定的，文獻[4]第24.4.3節有一些這方面的讨論，介紹了定性的trace plot方法和定量的Estimated potential scale reduction (EPSR)方法，可以做參考。
另外，在使用MCMC時，究竟是選擇一條很長的MC進行采樣以減少burn in時間還是選擇多條短的MC(優點？)，文獻[4]也給出了一個經驗選擇：用中等數量的中等長度的MC(比如3條MC，每條長度為\(10^{6}\))。
參考文獻
[1] Monte Carlo Strategies in Scientific Computing
[2] Numerical recipes
[3] Monte Carlo Statistical Methods
[4] Machine learning: a probabilistic perspective
![python運作mcmc為何老出錯_Python可視化解析MCMC[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)