轉自：http://blog.csdn.net/xianlingmao/article/details/7768833
引入
我們會遇到很多問題無法用分析的方法來求得精確解，例如由於式子特別，真的解不出來。這時就需要找一種方法求其近似解，並且有手段能測量出這種解的近似程度 （比如漸進性，上下限什么的）
隨機模擬的基本思想
現在假設我們有一個矩形的區域R（大小已知），在這個區域中有一個不規則的區域M（即不能通過公式直接計算出來），現在要求取M的面積？ 怎么求？近似的方法很多，例如：把這個不規則的區域M划分為很多很多個小的規則區域，用這些規則區域的面積求和來近似M，另外一個近似的方法就是采樣的方法，我們抓一把黃豆，把它們均勻地鋪在矩形區域，如果我們知道黃豆的總個數S，那么只要我們數數位於不規則區域M中的黃豆個數S1，那么我們就可以求出M的面積：M=S1*R/S。
在機器學習或統計計算領域，我們常常遇到這樣一類問題：即如何求取一個定積分：\inf _a ^b f(x) dx， 如歸一化因子等。
如何來求解這類問題呢？當然如果定積分可以解析求出，直接求出即可，如果不能解析求出，只能求取近似解了，常用的近似方法是采用蒙特卡洛積分，即把上述式子改寫為：
\inf _a^b f(x)*g(x)/g(x) dx = \inf _a^b (1/g(x)) *f(x)*g(x) dx
那么把f(x)/g(x)作為一個函數，而把g(x)看做是[a,b]上的一個概率分布，抽取n個樣本之后，上述式子可以繼續寫為：{\sum _1^n [f(x_i)/g(x_i)]}/n，當n趨向無窮大的時候，根據大數定理，上述式子是和要求的定積分式子是相等的，因此可以用抽樣的方法來得到近似解。
通過上述兩個例子，我們大概能夠理解抽樣方法解決問題的基本思想，其基本思路就是要把待解決的問題轉化為一種可以通過某種采樣方法可以解決的問題，至於怎么轉化，還是挺有創造性，沒有定法。因此隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。
常見的抽樣方法
直接抽樣法
因為較為簡單，而且只能解決很簡單的問題，一般是一維分布的問題
接受-拒絕抽樣（Acceptance-Rejection sampling)
為了得到一個分布的樣本，我們通過某種機制得到了很多的初步樣本，然后其中一部分初步樣本會被作為有效的樣本（即要抽取的分布的樣本），一部分初步樣本會被認為是無效樣本舍棄掉。這個算法的基本思想是：我們需要對一個分布f(x)進行采樣，但是卻很難直接進行采樣，所以我們想通過另外一個容易采樣的分布g(x)的樣本，用某種機制去除掉一些樣本，從而使得剩下的樣本就是來自與所求分布f(x)的樣本。
它有幾個條件：1）對於任何一個x，有f(x)<=M*g(x); 2) g(x)容易采樣；3) g(x)最好在形狀上比較接近f(x)。具體的采樣過程如下：
1. 對於g(x)進行采樣得到一個樣本xi, xi ~ g(x);
2. 對於均勻分布采樣 ui ~ U(a,b);
3. 如果ui<= f(x)/[M*g(x)], 那么認為xi是有效的樣本；否則舍棄該樣本； （# 這個步驟充分體現了這個方法的名字：接受-拒絕）
4. 反復重復步驟1~3，直到所需樣本達到要求為止。
這個方法可以如圖所示：
(說明：這是從其他地方弄來的圖，不是自己畫的，符號有些和文中不一致，其中\pi(x) 就是文中的f(x)，q(x)就是文中的g(x) )
重要性抽樣(Importance sampling)
重要性采樣和蒙特卡洛積分密切相關，看積分：
\inf f(x) dx = \inf f(x)*(1/g(x))*g(x) dx, 如果g(x)是一個概率分布，從g(x)中抽取N個樣本，上述的式子就約等於(1/N)* \sum f(xi)*(1/g(xi))。這相當於給每個樣本賦予了一個權重，g(xi)大意味着概率大，那么N里面含有這樣的樣本xi就多，即這些樣本的權重大，所以稱為重要性抽樣。
抽樣的步驟如下：
1. 選擇一個容易抽樣的分布g(x), 從g(x)中抽取N個樣本；
2. 計算(1/N)* \sum f(xi)*(1/g(xi))，從而得到近似解。（pku，sewm，shinning）
MCMC抽樣方法
無論是拒絕抽樣還是重要性采樣，都是屬於獨立采樣，即樣本與樣本之間是獨立無關的，這樣的采樣效率比較低，如拒絕采樣，所抽取的樣本中有很大部分是無效的，這樣效率就比較低，MCMC方法是關聯采樣，即下一個樣本與這個樣本有關系，從而使得采樣效率高。MCMC方法的基本思想是：通過構建一個markov chain使得該markov chain的穩定分布是我們所要采樣的分布f(x)。如果這個markov chain達到穩定狀態，那么來自這個chain的每個樣本都是f(x)的樣本，從而實現抽樣的目的。
Metropolis-Hasting算法
假設要采樣的概率分布是\pi(x)，現在假設有一個概率分布p(y|x)，使得\pi(x)*p(y|x) = \pi(y)*p(x|y)成立，稱細致平衡公式，這個細致平衡公式是markov chain能達到穩定分布的必要條件。因此關鍵是構建出一個概率分布p(y|x)使得它滿足細致平衡。現在假設我們有一個容易采樣的分布q(y|x)（稱為建議分布)，對於目前的樣本x，它能夠通過q(y|x)得到下一個建議樣本y，這個建議樣本y按照一定的概率被接受或者不被接受，稱為比率\alpha(x, y) = min{1, q(x|y)*\pi(y)/[q(y|x)*\pi(x)]}。即如果知道樣本xi，如何知道下一個樣本x_{i+1}是什么呢？就是通過q(y|xi)得到一個建議樣本y，然后根據\alpha(xi, y)決定x_{i+1}=y 還是x_{i+1}=xi。可以證明分布q(y|x)*\alpha(x,y)滿足細致平衡，同時可以證明這樣抽取得到的樣本是分布\pi(x)的樣本。具體的步驟如下：
1. 給定一個起始樣本x_0和一個建議分布q(y|x)；
2. 對於第i個樣本xi，通過q(y|xi)得到一個建議樣本y；計算比率\alpha(xi, y)= min{1, q(xi|y)*\pi(y)/[q(y|xi)*\pi(xi)]};
3. 抽取一個均勻分布樣本ui ~ U(0,1)，如果ui <= \alpha(xi,y)，則x_{i+1} = y；否則x_{i+1} = xi；
4. 重復步驟2~3，直到抽取到想要的樣本數量為止。
如果，建議分布q(y|x) 滿足：q(y|x) = q(x|y)，即對稱，這個時候比率\alpha(x, y) = min{1, \pi(y)/\pi(x)}就是1953年最原始的算法，后來hasting把這個算法擴展了，不要求建議分布式對稱的，從而得到了上述的算法。然而這個算法有一個缺點，就是抽樣的效率不高，有些樣本會被舍棄掉。從而產生了Gibbs算法。
Gibbs采樣算法
Gibbs算法，很簡單，就是用條件分布的抽樣來替代全概率分布的抽樣。例如，X={x1,x2,...xn}滿足分布p(X)，如何對p(X)進行抽樣呢？如果我們知道它的條件分布p(x1|X_{-1}),...,p(xi|X_{-i}),....，其中X_{-i}表示除了xi之外X的所有變量。如果這些條件分布都是很容易抽樣的，那么我們就可以通過對條件分布的抽樣來對全概率分布p(X)進行抽樣。
Gibbs采樣算法的步驟：
1. 給定一個初始樣本X0={x10,x20,...,xn0}
2.已知一個樣本Xi={x1i,x2i,...,xni}，對於x1_{i+1}進行抽樣，x1_{i+1} ~ p(x1|Xi_{-1})
3. 對於x2_{i+1}進行抽樣，x2_{i+1} ~ p(x2|x1_{i+1}, x3i,...xni)
................................................................
4.對於xn_{i+1}進行抽樣，xn_{i+1} ~ p(xn|x1_{i+1}, x2_{i+1},...x_{n-1}_{i+1})
5.步驟2~4可以得到X的一個樣本，然后重復步驟2~4可以不斷地得到X的樣本。
當然無論是metropolis-hasting算法還是gibbs算法，都有一個burn in的過程，所謂burn in的過程就是因為這個兩個算法本身都是markov chain的算法，要達到穩定狀態需要一定的步驟才能達到，所以需要一個burn in過程，只有在達到平衡狀態時候得到的樣本才能是平衡狀態時候的目標分布的樣本，因此，在burn in過程中產生的樣本都需要被舍棄。如何判斷一個過程是否達到了平衡狀態還沒有一個成熟的方法來解決，目前常見的方法是看是否狀態已經平穩（例如畫一個圖，如果在較長的過程中，變化已經不大，說明很有可能已經平衡）當然這個方法並不能肯定一個狀態是否平衡，你可以舉出反例，但是卻是實際中沒有辦法的辦法。