1.隨機模擬

隨機模擬，又名蒙特卡洛方法（Monte Carlo Simulation）,發展始于20世紀40年代。
現代的統計模擬方法最早由數學家烏拉姆提出，被Metropolis命名為蒙特卡洛方法。蒙特卡洛是著名的賭場，賭博總是和統計密切相關的。布豐當年用于計算n的著名的投針試驗就是蒙特卡羅模擬實驗，隨機采樣的方法其實數學家們很早就知道，但是在計算機出現以前，隨機數生成的成本很高，所以該方法也沒有實用價值，隨著計算機技術的迅猛發展，隨機模擬技術很快進入實用階段，對那些用確定算法不可行或不可能解決的問題，蒙特卡洛方法常常為人們帶來希望。
統計模擬中有一個重要的問題就是給定一個概率分布p(x)，我們如何在計算機中生成它的樣本。一般而言，均勻分布Unifrom(0,1)的樣本是相對容易生成的。通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成[0，1]之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布Uniform(0,1)的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。

其它常見的概率分布，無論是連續的還是離散的分布，都給予Unifrom（0，1）樣本生成，正態分布可以通過著名的Box-Muller變換得到。
【正態分布】若隨機變量X服從一個數學期望為μ、方差為σ^2的正態分布，記為N(μ，σ^2)。其概率密度函數為正態分布的期望值μ決定了其位置，其標準差σ決定了分布的幅度。當μ = 0,σ = 1時的正態分布是標準正態分布。
正態曲線呈鐘型，兩頭低，中間高，左右對稱因其曲線呈鐘形，因此人們又經常稱之為鐘形曲線。

【Box-Muller變換】
如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，

則 Z0,Z1 獨立且服從標準正態分布
正態值 Z 有一個等于 0 的平均值和一個等于 1 的標準偏差，可使用以下等式將 Z 映射到一個平均值為 m、標準偏差為 sd 的統計量 X：X = m + (Z * sd)

其它幾個著名的連續分布，包括指數分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通過類似的數學變換得到；離散的分布通過均勻分布更加容易生成。更多的統計分布如何通過均勻分布的變換生成出來，大家可以參考統計計算的書，其中 Sheldon M. Ross 的《統計模擬》是寫得非常通俗易懂的一本。
不過我們并不是總是這么幸運的，當p(x)的形式很復雜，或者 p(x) 是個高維的分布的時候，樣本的生成就可能很困難了。 譬如有如下的情況

p(x)=p~(x)∫p~(x)dx,而 p~(x) 我們是可以計算的，但是底下的積分式無法顯式計算。

p(x,y) 是一個二維的分布函數，這個函數本身計算很困難，但是條件分布 p(x|y),p(y|x)的計算相對簡單;如果 p(x) 是高維的，這種情形就更加明顯。

此時就需要使用一些更加復雜的隨機模擬的方法來生成樣本。而本節中將要重點介紹的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一種，這兩個方法在現代貝葉斯分析中被廣泛使用。要了解這兩個算法，我們首先要對馬氏鏈的平穩分布的性質有基本的認識。

1.2蒙特卡洛數值積分

如果我們要求f(x)的積分，如

而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：

這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：

因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。

1.3 Monte Carlo principle

Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。

1.4 接受-拒絕抽樣（Acceptance-Rejection sampling)

很多實際問題中，p(x)是很難直接采樣的的，因此，我們需要求助其他的手段來采樣。既然 p(x) 太復雜在程序中沒法直接采樣，那么我設定一個程序可抽樣的分布 q(x) 比如高斯分布，然后按照一定的方法拒絕某些樣本，達到接近 p(x) 分布的目的，其中q(x)叫做 proposal distribution 。

具體操作如下，設定一個方便抽樣的函數 q(x)，以及一個常量 k，使得 p(x) 總在 kq(x) 的下方。（參考上圖）

• x 軸方向：從 q(x) 分布抽樣得到 a。(如果是高斯，就用之前說過的 tricky and faster 的算法更快）
• y 軸方向：從均勻分布（0, kq(a)) 中抽樣得到 u。
• 如果剛好落到灰色區域： u > p(a), 拒絕， 否則接受這次抽樣
• 重復以上過程

在高維的情況下，Rejection Sampling 會出現兩個問題，第一是合適的 q 分布比較難以找到，第二是很難確定一個合理的 k 值。這兩個問題會導致拒絕率很高，無用計算增加。

1.5 重要性抽樣(Importance sampling)

Importance Sampling 也是借助了容易抽樣的分布 q (proposal distribution)來解決這個問題，直接從公式出發：

很多實際問題中，p(x)是很難直接采樣的的，因此，我們需要求助其他的手段來采樣。既然 p(x) 太復雜在程序中沒法直接采樣，那么我設定一個程序可抽樣的分布 q(x) 比如高斯分布，然后按照一定的方法拒絕某些樣本，達到接近 p(x) 分布的目的，其中q(x)叫做 proposal distribution 。

第一個圖表示 p 分布， 第二個圖的陰影區域 f = 1，非陰影區域 f = 0, 那么一個良好的 q 分布應該在左邊箭頭所指的區域有很高的分布概率，因為在其他區域的采樣計算實際上都是無效的。這表明 Importance Sampling 有可能比用原來的 p 分布抽樣更加有效。

但是可惜的是，在高維空間里找到一個這樣合適的 q 非常難。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出現，要找到一個同時滿足 easy to sample 并且 good approximations 的 proposal distribution, it is often impossible！

1.6 馬氏鏈及其分布

氏鏈是指考察一個隨機過程，若己知現在t的狀態X(t)，那么將來的狀態X(t+n)取值（或取某些狀態）的概率與過去狀態X(s)（s小于t）取值無關，或更簡單的說，己知現在，將來與過去無關（條件獨立），則稱此性質為馬爾可夫性（無后效性或簡稱馬氏性）。

馬氏鏈的數學定義很簡單

也就是狀態轉移的概率只依賴于前一個狀態。
如圖：

【轉移概率矩陣】：矩陣各元素都是非負的，并且各行元素之和等于1，各元素用概率表示，在一定條件下是互相轉移的，故稱為轉移概率矩陣。如用于市場決策時，矩陣中的元素是市場或顧客的保留、獲得或失去的概率。P(k)表示k步轉移概率矩陣。

舉一個例子，如果當前狀態為 u(x) = (0.5, 0.2, 0.3), 那么下一個矩陣的狀態就是 u(x)T = (0.18, 0.64, 0.18), 依照這個轉換矩陣一直轉換下去，最后的系統就趨近于一個穩定狀態 (0.22, 0.41, 0.37) (此處只保留了兩位有效數字)。而事實證明無論你從那個點出發，經過很長的 Markov Chain 之后都會匯集到這一點。[2]

再舉一個例子，社會學家經常把人按其經濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1,2,3 分別代表這三個階層。社會學家們發現決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是 0.65, 屬于中層收入的概率是 0.28, 屬于上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉移概率如下


使用矩陣的表示方式，轉移概率矩陣記為



我們發現從第7代人開始，這個分布就穩定不變了，事實上，在這個問題中，從任意初始概率分布開始都會收斂到這個上面這個穩定的結果。

注：要求圖是聯通的（沒有孤立點），同時不存在一個聯通的子圖是沒有對外的出邊的（就像黑洞一樣）。

這個馬氏鏈的收斂定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以這個定理作為理論基礎的。

對于給定的概率分布p(x),我們希望能有便捷的方式生成它對應的樣本。由于馬氏鏈能收斂到平穩分布， 于是一個很的漂亮想法是：如果我們能構造一個轉移矩陣為P的馬氏鏈，使得該馬氏鏈的平穩分布恰好是p(x), 那么我們從任何一個初始狀態x0出發沿著馬氏鏈轉移, 得到一個轉移序列 x0,x1,x2,?xn,xn+1?,， 如果馬氏鏈在第n步已經收斂了，于是我們就得到了 π(x) 的樣本xn,xn+1?。

這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統的平穩性質， Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基于馬氏鏈的蒙特卡羅方法，即Metropolis算法，并在最早的計算機上編程實現。Metropolis 算法是首個普適的采樣方法，并啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis的這篇論文被收錄在《統計學中的重大突破》中， Metropolis算法也被遴選為二十世紀的十個最重要的算法之一。

我們接下來介紹的MCMC 算法是 Metropolis 算法的一個改進變種，即常用的 Metropolis-Hastings 算法。由上一節的例子和定理我們看到了，馬氏鏈的收斂性質主要由轉移矩陣P 決定, 所以基于馬氏鏈做采樣的關鍵問題是如何構造轉移矩陣P,使得平穩分布恰好是我們要的分布p(x)。如何能做到這一點呢？我們主要使用如下的定理。


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

1.7 MCMC——Gibbs Sampling算法

參考資料
[1]http://blog.csdn.net/xianlingmao/article/details/7768833

[2] http://www.cnblogs.com/daniel-D/p/3388724.html

[3] http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

[4] An Introduction to MCMC for Machine Learning，2003

[5] Introduction to Monte Carlo Methods

[6]http://www.cnblogs.com/xbinworld/p/4266146.html

2.隨機模擬的R語言實現

2.1Markov Chain （馬爾科夫鏈）

假設 f(t) 是一個時間序列，Markov Chain是假設f(t+1)只與f(t)有關的隨機過程。

Implement in R:

N = 10000 signal = vector(length = N) signal[1] = 0 for (i in 2:N) { # random select one offset (from [-1,1]) to signal[i-1] signal[i] = signal[i-1] + sample(c(-1,1),1) } plot( signal,type = 'l',col = 'red')

這里的問題：
在用RStudio運行時，會提示Error in plot.new() : figure margins too large。這個錯誤本質上來講就是畫的圖在畫布上展不開。這里有兩種解決方法：
1。拖動Rstudio的畫布，讓畫布的區域大一點。（簡單粗暴。。）
2。問題在于畫圖時的邊界空白行數（或列數）較大，因此報錯。
可以用par(“mar”)查看默認的邊界設置：
> par("mar")
[1] 5.1 4.1 4.1 2.1
然后在畫圖時，對“mar”參數進行用戶自定義就可以了

opar <- par(no.readonly = TRUE) par(mar = rep(2, 4)) # to be simple # Plot Code par(opar)

2.2Random Walk（隨機游走）

如布朗運動，只是上面Markov Chain的二維拓展版：
Implement in R:

N = 100 x = vector(length = N) y = vector(length = N) x[1] = 0 y[1] = 0 for (i in 2:N) { x[i] = x[i-1] + rnorm(1) y[i] = y[i-1] + rnorm(1) } plot(x,y,type = 'l', col='red')

2.3MCMC具體方法：

MCMC方法最早由Metropolis（1954）給出，后來Metropolis的算法由Hastings改進，合稱為M-H算法。M-H算法是MCMC的基礎方法。由M-H算法演化出了許多新的抽樣方法，包括目前在MCMC中最常用的Gibbs抽樣也可以看做M-H算法的一個特例[2]。

概括起來，MCMC基于這樣的理論，在滿足【平衡方程】（detailed balance equation）條件下，MCMC可以通過很長的狀態轉移到達穩態。
【平衡方程】：
pi(x) * P(y|x) = pi(y) * P(x|y)

其中pi指分布，P指概率。這個平衡方程也就是表示條件概率（轉化概率）與分布乘積的均衡.

3.1 M-H法
1. 構造目標分布，初始化x0
2. 在第n步，從q(y|x_n) 生成新狀態y
3. 以一定概率（(pi(y) * P(x_n|y)) / (pi(x) * P(y|x_n))）接受y < PS: 看看上面的平衡方程，這個概率表示什么呢？參考這里和[1]>

implementation in R:

N = 10000 x = vector(length = N) x[1] = 0 # uniform variable: u u = runif(N) m_sd = 5 freedom = 5 for (i in 2:N) { y = rnorm(1,mean = x[i-1],sd = m_sd) print(y) #y = rt(1,df = freedom) p_accept = dnorm(x[i-1],mean = y,sd = abs(2*y+1)) / dnorm(y, mean = x[i-1],sd = abs(2*x[i-1]+1)) #print (p_accept) if ((u[i] <= p_accept)) { x[i] = y print("accept") } else { x[i] = x[i-1] print("reject") } } plot(x,type = 'l') dev.new() hist(x)

(這里和參考代碼相比少了一幅圖，還在看為什么。。。。)

3.4Gibbs采樣

第n次，Draw from ，迭代采樣結果接近真實p

也就是每一次都是固定其他參數，對一個參數進行采樣。比如對于二元正態分布，其兩個分量的一元條件分布仍滿足正態分布：

那么在Gibbs采樣中對其迭代采樣的過程，實現如下：

#define Gauss Posterior Distribution p_ygivenx <- function(x,m1,m2,s1,s2) { return (rnorm(1,m2+rho*s2/s1*(x-m1),sqrt(1-rho^2)*s2 )) } p_xgiveny <- function(y,m1,m2,s1,s2) { return (rnorm(1,m1+rho*s1/s2*(y-m2),sqrt(1-rho^2)*s1 )) } N = 5000 K = 20 #iteration in each sampling x_res = vector(length = N) y_res = vector(length = N) m1 = 10; m2 = -5; s1 = 5; s2 = 2 rho = 0.5 y = m2 for (i in 1:N) { for(j in 1:K) { x = p_xgiveny(y, m1,m2,s1,s2) y = p_ygivenx(x, m1,m2,s1,s2) # print(x) x_res[i] = x; y_res[i] = y; } } hist(x_res,freq = 1) dev.new() plot(x_res,y_res) library(MASS) valid_range = seq(from = N/2, to = N, by = 1) MVN.kdensity <- kde2d(x_res[valid_range], y_res[valid_range], h = 10) #估計核密度 plot(x_res[valid_range], y_res[valid_range], col = "blue", xlab = "x", ylab = "y") contour(MVN.kdensity, add = TRUE)#二元正態分布等高線圖 #real distribution real = mvrnorm(N,c(m1,m2),diag(c(s1,s2))) dev.new() plot(real[1:N,1],real[1:N,2])

Reference:
1.http://www2.isye.gatech.edu/~brani/isyebayes/bank/handout10.pdf
2.http://site.douban.com/182577/widget/notes/10567181/note/292072927/
3. book: http://statweb.stanford.edu/~owen/mc/
4.Classic：http://cis.temple.edu/~latecki/Courses/RobotFall07/PapersFall07/andrieu03introduction.pdf
5.https://segmentfault.com/q/1010000000575941
6.http://blog.csdn.net/abcjennifer/article/details/25908495

3 關于隨機模擬和MCMC算法的個人總結和感悟

MCMC算法作為二十世紀以來最好的十大算法之一，有其突出的獨特性和精妙性，被用在很多實際應用當中。

在看這個算法之前，我對統計學的相關知識和算法并沒有全面和系統的了解，很多知識都是模棱兩可大而化之的。所以這也導致我在算法的理論部分花費了很多時間，從一開始簡單的針對性的去看這個算法本身，到為了更好的理解算法的思路而擴展開來了解更多的相關知識，并借此認真的把很多基礎的概念重新整理和重溫了一遍，雖然現在還是還有很多不甚明白的地方，所以以后還是要在這個領域繼續深入長期的學習。

在代碼實現部分，因為之前沒有涉及過R語言，獨立上手有些困難，所以代碼主要參考了網上的部分代碼，自己進行了再次實現，并對代碼句段進行了個人理解。在這個過程中，也總結了一些作為一個新手在代碼實現中的小問題，大部分也通過查資料和咨詢同學得到了解決。目前自己也已經準備從頭開始，認真細致的開始學習R語言。

在查詢資料的過程中，我發現網上關于MCMC的資料并不多，而且內容很零散，針對性也比較弱，學習時會花費大量的時間。所以我站在一個有一定數學基礎但是完全新手的角度上，對我能找到的一些資料（包括理論和實踐部分）進行了再處理和總結，并按照我的理解對文章內容進行了重新整合和糾錯，補充了相關基礎知識和注意要點，最終形成了這篇學習報告。文章中提到的內容所涉及到的參考資料，也都在各部分的最后附上了鏈接。

文章內容也會隨著我將來的深入學習，不斷的進行更新和再整理。

總結

以上是生活随笔為你收集整理的MCMC算法学习总结的全部內容，希望文章能夠幫你解決所遇到的問題。

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