一、MCMC 簡介
1. Monte Carlo 蒙特卡洛
蒙特卡洛方法（Monte Carlo）是一種通過特定分佈下的隨機數（或偽隨機數）進行模擬的方法。典型的例子有蒲豐投針、定積分計算等等，其基礎是大數定律。
蒙特卡洛方法有哪些優缺點如下：
- 優點：計算準確性由取樣的均勻程度決定；大大簡化問題複雜性
- 缺點：
- 由於要進行大量的抽樣計算，對計算機速度依賴性強
- 目前絕大多數隨機數發生器均為偽隨機數，一定程度上有偏
- 定積分求解問題中，對於\(\color{blue}{複雜或者高維的分佈}\)，利用蒙特卡洛方法生成隨機樣本非常困難
蒙特卡羅方法在金融工程學，巨集觀經濟學，生物醫學，計算物理學(如粒子輸運計算、量子熱力學計算、空氣動力學計算、核工程)等領域應用廣泛。
2. 隨機數生成演算法簡介
2.1 均勻分佈隨機數生成
均勻分佈隨機數生成最普遍的演算法為線性同餘演算法(LCG)，他是一個偽隨機數發生器，預設的初始種子為系統時間戳，演算法遞推公式如下：
\[
X_{i+1} = (A \times X_{i} + B) mod M
\]
其中：
- \(B, M\)互質
- M的所有質因子的積能整除\((A-1)\)：比如M是3的倍數則\((A-1)\)也是
- \(A, B, X_0\)都比\(M\)小
- \(A,B\)是正整數
python程式碼如下:
import matplotlib.pyplot as plt
import time
import math
def random(num = 2000):
M = 1000
A = 3^4
B = 1
x_0 = time.time()
result = [x_0]
while len(result) < num:
result.append((result[-1]*A+B) % M)
result.pop(0)
return result
result = random(2000)
xaixs = [i for i,j in enumerate(result)]
print result[:10]
f = plt.figure(1)
p1 = plt.scatter(xaixs,result)
#p2 = plt.hist(result)
plt.show()
目前純的隨機數生成演算法非常的複雜，而且效能消耗非常的大，所以在此不做介紹。
2.2 標準正態分佈隨機數生成
已知均勻分佈隨機數如何生成，我麼可以根據其生成標準正態分佈隨機數，Box和Muller(1958) (G. E. P. Box and Mervin E. Muller, A Note on the Generation of Random Normal Deviates, The Annals of Mathematical Statistics (1958), Vol. 29, No. 2 pp. 610–611)提出了產生\(N(0,1)\)隨機數的簡單方法：設\(U_1,U_2\)是獨立同分布的\(U(0,1)\)變數，記:
\[ \begin{cases} X_1 = (-2 lnU_1)^{1/2}cos(2\pi U_2) \\ X_2 = (-2 lnU_2)^{1/2}cos(2\pi U_1) \end{cases} \]
則\(X_1\)和\(X_2\)相互獨立且均服從標準正態分佈。python程式碼如下：
import matplotlib.pyplot as plt
import time
import math
def rnorm(num = 1000):
uniform1 = random(num)
uniform2 = random(num)
return [math.sqrt(-2*math.log(i/1000,math.e))*math.cos(2*math.pi*j/1000) for i,j in zip(uniform1,uniform2)]
result = rnorm()
xaixs = [i for i,j in enumerate(result)]
f = plt.figure(1)
#p1 = plt.scatter(xaixs,result)
p2 = plt.hist(result)
plt.show()2.3 其他簡單分佈隨機數生成
對於已知pdf或是cdf的分佈，可以根據cdf求反函式的方法來生成該分佈的隨機數。因為cdf值域為\([0,1]\)，我們只需要獲得均勻分佈隨機數\(U_1 \sim U(0,1)\)，即可根據反函式\(cdf^{-1}\)得到此次抽得的隨機數。對於再複雜一點的分佈，我們可以通過構造包絡分佈進行Rejection Sampling，如ARS和ARMS演算法等等。
3. Markov Chain 馬爾科夫鏈
在蒙特卡洛模擬中，我們在後驗分佈中抽取樣本，當這些樣本獨立時，利用大數定律樣本均值會收斂到期望值。如果得到的樣本是不獨立的，那麼就要藉助於馬爾科夫鏈進行抽樣。馬爾科夫鏈又稱為馬爾科夫過程是一種離散的隨機過程，隨機過程可以看做一個隨時間變化的隨機變數序列。
對一個馬爾科夫鏈來說，未來狀態只與當前t時刻有關，而與t時刻之前的歷史狀態無關。馬爾科夫鏈的每個狀態的狀態轉移概率構成了其狀態轉移矩陣\(P\)
** 當步長\(n\)足夠大時，一個非週期且任意狀態聯通的馬氏鏈可以收斂到一個平穩分佈\(\pi(x)\) **，馬氏鏈收斂定理是MCMC的基礎。
二、Sampling Algorithm 抽樣演算法
1. MCMC與M-H演算法
MCMC方法就是構造合適的馬爾科夫鏈進行抽樣而使用蒙特卡洛方法進行積分計算。 既然馬爾科夫鏈可以收斂到平穩分佈。我們可以建立一個以\(\pi\)為平穩分佈的馬爾科夫鏈，對這個鏈執行足夠長時間之後，可以達到平穩狀態。此時馬爾科夫鏈的值就相當於在分佈\(\pi(x)\)中抽取樣本。利用馬爾科夫鏈進行隨機模擬的方法就是MCMC。
MCMC的主要問題是：如何構造轉移矩陣\(P\),使得平穩分佈恰好是我們要的分佈\(\pi (x)\)。由此，我們首先有了M-H演算法，它用到可馬爾科夫鏈的另一個性質，如果具有轉移矩陣\(P\)和分佈\(\pi(x)\)的馬氏鏈對所有的狀態\(i,j\)滿足下面的等式：
\[
\pi(i)P(i,j) = \pi(j)P(j,i)
\]
這個等式稱為細緻平衡方程，此時該馬氏鏈的分佈\(\pi(x)\)是平穩的。一般情況下，任意給的\(\pi(x)\)不一定平穩，那麼可以建構函式\(\alpha(i,j)\),使得：
\[ \pi(i)P(i,j)\alpha(i,j) = \pi(j)P(j,i)\alpha(j,i) \]
此時取\(\alpha(i,j) = \pi(j)P(j,i)\)，\(\alpha(j,i) = \pi(i)P(i,j)\)時等式(2)必然成立。
(2)式中的\(\alpha(i,j)\)又叫接受率，是接受狀態從\(i\)轉移到\(j\)的概率。如果\(\alpha\)過小，將會導致頻繁拒絕轉移，馬氏鏈難以收斂，所以我們考慮放大\(\alpha\)，使得\(\alpha(i,j)\)和\(\alpha(j,i)\)中大的值為\(1\)，小的值同比放大為\(\frac{\pi(j)P(j,i)}{\pi(i)P(i,j)}\)，那麼\(\alpha(i,j) = min\{1,\frac{\pi(j)P(j,i)}{\pi(i)P(i,j)}\}\)。在馬爾科夫鏈狀態轉移過程中，我們以概率\(\alpha(i,j)\)接受轉移，從而得到最終的穩定分佈，這種演算法稱為M-H演算法。
M-H演算法需要提前構造提議分佈(proposal distribution)\(\pi(x)\)，然後經過提議分佈取樣後計算出接受概率，對取樣進行接受或是拒絕，最後得到服從特定分佈的樣本。優點是通過構造精妙的馬爾科夫鏈，提高了抽樣的效率，尤其是面對條件分佈抽樣時。缺點是需要找到合適的提議分佈，而且收斂速度可能較慢。
現欲對自由度為5的t分佈抽取10000個樣本，我們分別用M-H演算法和R語言自帶的t分佈抽樣命令rt進行抽樣，然後對比二者結果。
對於M-H抽樣，我們選擇標準正態分佈為其提議分佈，進行10000次模擬，R語言程式碼和得到結果如下：
ptm <- proc.time()
#### M-H algorithm ####
N = 10000
x = vector(length = N)
x[1] = 0
# uniform variable: u
u = runif(N)
sd = 1
df = 5
# sample the student distrubution with 5 freedom
for (i in 2:N) {
y = rnorm(1,mean = x[i-1],sd = sd)
p_accept = dt(y,df) * pnorm(x[i-1],y,sd) / pnorm(x[i-1],y,sd) / (dt(x[i-1],df))
if ((u[i] <= p_accept))
{ x[i] = y }
else
{ x[i] = x[i-1] }
}
HM_result <- x
proc.time()-ptm
#plot(x,type = 'l')
test = rt(10000,5)
par(mfrow=c(2,2))
hist(test)
plot(density(test),main = "density of test dataset")
hist(HM_result)
plot(density(HM_result),main = "density of MH result") 
2. Gibbs抽樣
M-H抽樣在較低維度值運算中是比較方便的，然而在維數較高時，選擇適當的提議分佈並非易事，而Gibbs抽樣僅僅需要了解滿條件分佈(full conditi- onal distribution)，因此在高維積分時更具有優勢（朱慧明 等, 《貝葉斯計量經濟模型》, 2009,科學出版社）。
n維變數的Gibbs抽樣方法如下:給定抽樣初始點\(x^{(0)} = (x_1^{(0)},x_2^{(0)},...,x_n^{(0)}\) \()\)，並假設第\(t-1\)次迭代的迭代值為\(x^{(t-1)}\)，則第\(t\)次迭代由以下步驟完成：
- 由滿條件分佈\(\pi(x_1 | x_2^{(t-1)},...,x_n^{(t-1)})\)抽取\(x_1^{(t)}\)
- 由滿條件分佈\(\pi(x_2 | x_1^{(t-1)},x_3^{(t-1)}...,x_n^{(t-1)})\)抽取\(x_2^{(t)}\)
- ... ...
- 由滿條件分佈\(\pi(x_n | x_1^{(t-1)},...,x_{n-1}^{(t-1)})\)抽取\(x_n^{(t)}\)
至此完成\(t-1\)到\(t\)的一次迭代，重複迭代直到\(x^{(t)}\)分佈穩定後即可。
Gibbs抽樣例項
以二維正態分佈抽樣為例，不妨設\(X = (X_1,X_2)\)服從二元正態分佈，即:
\[
\begin{eqnarray}
x \sim N
\left(
\begin{pmatrix}
0 \\ 0
\end{pmatrix}
,
\begin{pmatrix}
1 & \rho \\
\rho & 1
\end{pmatrix}
\right),\ 0<\rho\le1
\end{eqnarray}
\]
則有:
- $x_2^{(t+1)}| x_1^{(t+1)} \sim N(\rho x_1^{t+1},1-\rho^2) $
- \(x_1^{(t+1)} | x_2^{(t+1)} \sim N(\rho x_2^{t+1},1-\rho^2)\)
並且可以證明：
\[
\begin{eqnarray}
\begin{pmatrix}
x_1^{(t)} \\
x_2^{(t)}
\end{pmatrix}
\sim
N \left(
\begin{pmatrix}
\rho^{2t-1}x_1^{(0)} \\
\rho^{2t}x_2^{(0)}
\end{pmatrix}
,
\begin{pmatrix}
1 - \rho^{4t-2} & \rho - \rho^{4t-1} \\
\rho - \rho^{4t-1} & 1 - \rho^{4t}
\end{pmatrix}
\right)
\end{eqnarray}
\]
當\(t \to \infty\)時，\((x_1^{t},x_2^{t})\)的分佈將收斂到目標分佈。R語言程式碼和結果如下：
n <- 5000 #抽樣個數
n_1 <- 2500 # 收斂個數
X <- matrix(0, n, 2)
mu1 <- 0
mu2 <- 0
sigma1 <- 1
sigma2 <- 1
rho <- 0.5
s1.c <- sqrt(1 - rho^2) * sigma1
s2.c <- sqrt(1 - rho^2) * sigma2
X[1, ] <- c(mu1, mu2) #初始化
for (i in 2:n) {
x2 <- X[i - 1, 2]
m1.c <- mu1 + rho * (x2 - mu2) * sigma1/sigma2
X[i, 1] <- rnorm(1, m1.c, s1.c)
x1 <- X[i, 1]
m2.c <- mu2 + rho * (x1 - mu1) * sigma2/sigma1
X[i, 2] <- rnorm(1, m2.c, s2.c)
}
b <- n_1 + 1
x.mcmc <- X[b:n, ]
library(MASS)
MVN.kdensity <- kde2d(x.mcmc[, 1], x.mcmc[, 2], h = 5) #估計核密度
par(mfrow=c(2,2))
plot(x.mcmc, col = "blue", xlab = "X1", ylab = "X2") #繪製二元正態分佈等高線圖
contour(MVN.kdensity, add = TRUE)
hist(x.mcmc[, 1],main = "histogram of x1")
hist(x.mcmc[, 2],mian = "histogram of x2")
x<-y<-seq(-4,4,length=20)
f<-function(x,y){(exp(-0.5*x^2-0.5*y^2))/(2*pi)}
z<-outer(x,y,f)
persp(x,y,z,theta=45,phi=25,col='lightblue') 
三、MCMC和線性迴歸模型
當我們討論貝葉斯估計方法和線性迴歸模型時，會包括:
- \(\beta\)的先驗分佈
- \(\sigma\)的先驗分佈
給定一元線性迴歸模型:
\[
Y_i = \beta_0 + \beta_1 X_i + \epsilon_i
\]
我們將其中心化，變為：
\[
\begin{eqnarray}
Y_i &=& \beta_0^{*} + \beta_1 (X_i - \bar{X}) + \epsilon_i \\
Y_i &=& \beta_0^{*} - \beta_1 \bar{X} + \beta_1 X_i + \epsilon_i
\end{eqnarray}
\]
這樣變化後，我們發現\(\beta_1\)是不變的，由於\(X\)的擾動沒變而\(\bar{X}\)無擾動。而\(\beta_0\)的期望相對增加了，有\(E(\beta_0) = \bar{Y}\)，且自變數\(X_i - \bar{X}\)均值為0。這樣模型的極大似然值可以寫成：
\[
y_i | x_i,\beta_0,\beta_1,\sigma^2 \sim N(\beta_0 + \beta_1(x_i - \bar{x}),\sigma^2),i = 1,2,...,n
\]
貝葉斯模型充分利用樣本的先驗資訊，對於中心化一元線性模型(5)式， 我們有先驗資訊：
\[
p(\beta_0, \beta_1,\sigma^2) \propto \frac{1}{\sigma^2},-\infty < \beta_0,\beta_1 < \infty,0<\sigma^2<\infty
\]
同時我們有三個充分統計量：
\[
\begin{eqnarray}
\nonumber
\hat{\beta_0} &=& \bar{y} \\
\nonumber
\hat{\beta_1} &=& \frac{\sum_{i} (x_i - \bar{x})(y_i - \bar{y})}{\sum_i (x_i - \bar{x})^2} \\
\nonumber
s^2 & =& \frac{SSR}{n-2} = \frac{\sum_{i} (y_i - \hat{\beta_0} - \hat{\beta_1}(x_i - \bar{x})^2)}{n-2}
\end{eqnarray}
\]
根據這三個充分統計量、先驗分佈、極大似然函式，我們可以求得聯合後驗分佈\(p(\beta_0, \beta_1,\sigma^2|y)\)\footnote{Applied Bayesian Statistics with R and OpenBUGS, 187-190}
之後求得三個引數各自的邊緣後驗密度分佈:
\[
\begin{eqnarray}
\beta_0|y &\sim& t(\hat{\beta_0},\frac{s^2}{n},n-2) \\
\beta_1|y &\sim& t(\hat{\beta_1},\frac{s^2}{\sum (x_i - \bar{x})^2},n-2) \\
\sigma^2 | y &\sim& IG (\frac{n-2}{2},\frac{(n-2)s^2}{2})
\end{eqnarray}
\]
其中\(\beta_0,\beta_1\)的後驗分佈服從含位置和尺度引數的t分佈\footnote{if \(X \sim t(\mu,\sigma^2,freedom)\),then \((X - \mu) / \sigma^2 \sim t(freedom)\)}；而\(\sigma^2\)服從逆伽馬分佈。
得到引數的後驗分佈之後，我們就可以利用MCMC進行模擬並且收斂到最後的估計值。具體步驟如下：
- 引數賦予初始值(直接影響收斂速度) \(\beta_0^{(0)},\beta_1^{(0)},\sigma^2{(0)}\)
- 根據公式(5)更新\(\beta_0\)，其中\(\sigma^2 = (n-2) s /n\)
- 根據公式(6) 和公式(7)更新\(\beta_1,\sigma^2\)
- 重複(2) - (4) N輪，並觀察三個引數的收斂情況，從而得到穩健的估計值