采樣方法
實際應用中，經常需要獲得服從某一分布的樣本集。不過，手動生成一般來說不太現實，需要求助於計算機，而計算機則只能實現對均勻分布進行抽樣。其他的分布，甚至如高斯分布都是無法實現的。不過，通過均勻分布，可間接地生成服從其他分布的樣本。這點很重要，下面會看到，所有的隨機模擬都從均勻分布開始的，然后經過所需分布的約束，獲得所需樣本的。
Inverse CDF
最簡單，最直觀的方法是Inverse CDF，也稱為Inverse transform sampling。 借助於PDF（概率密度函數）與CDF（累積分布函數）的關系進行抽樣。
舉個簡單例子，設某一概率密度函數為：
這個密度函數的累積分布函數可以很容易的求得：
假設有一個樣本集來自此分布函數，那么，將每個樣本帶入累積分布函數，都可以得到對應的\(F(x)\)的值。如果樣本集來自隨機采樣，與其對應的累積分布函數值也會是隨機的，且服從(0,1)上的均勻分布。於是，為了得到這個的樣本，可以將這個過程顛倒過來，即在(0,1）的均勻分布上隨機采樣，然后帶入累積分布函數的反函數中，即可得到服從此分布的樣本 。
比如，\(F(x) = 0.8\)，我們需要求出其中的x是多少，這個x即是服從上述分布的樣本。求取的思路也很直觀，先求出\(F(x)\)的反函數，然后將0.8帶入，即可求出x的值是多少。上述累積分布反函數：
當 u =0.8時，帶入可得 \(x = F'(u) \approx 1.386\). 這個x即是服從上述分布的樣本，重復這個過程即可得到所需樣本集（如圖1 d圖）。
圖1 （a)概率密度函數；(b) 累積分布函數；(c) 累積分布函數的反函數；(d) 樣本集的分布。
可見，此種方法適用於可顯式地求得累積分布函數的反函數的分布。一般地，Inverse transform sampling 方法的抽樣過程可歸納為：
Given Cumulative Distribution Function of interest cdf(u), get its inverse function: icdf(x)
samples = list()
while True:
sample x from uniform(0,1)
y = icdf(x)
samples.append(y)
if satisfy some sampling condition:
break
對於常見的分布，如高斯分布不易求得其累積分布函數的常見分布，前人也給出了解決方法，如Box-Muller變換。
但對於不常見的分布，Inverse transform sampling應用起來就不太方便，因為累積分布函數的反函數不易求，甚至無法顯式地求出。這時就需要用其他方法，比如下面提到的接受-拒絕采樣（Acceptance-Rejection Sampling), 重要性采樣等等.
接受-拒絕采樣(Acceptance-Rejection Sampling)
當對某一分布\(p(x)\)直接抽樣比較困難時， 可以通過對另一相對容易的分布\(q(x)\)進行抽樣，然后保留其中服從\(p(x)\)的樣本，而剔除不服從\(p(x)\)的無效樣本。\(q(x)\)稱為建議分布函數（proposal distribution）。
例如，要對beta(2, 5)進行抽樣，beta分布的反函數不易求得，上述的反函數方法不太適用。那么可以引入一個相對容易抽樣的分布，如正態分布，甚至均勻分布等，這里不妨使用正態分布，比如N(0.2,1)。如上所述，抽樣可轉化為平面上的積分，需要proposal distribution 包裹或者說‘罩住’待抽樣分布，為此需將proposal distribution乘以一個系數 k，這里取k=2，具體如下圖：
接下來就可以進行采樣了：
-
從q(x) （proposal distribution，這里為N(0.2,0.3)隨機地抽取一個樣本, \(x_i\);
-
從均勻分布U(0,1)中隨機抽取一個樣本 \(u_i\);
-
比較\(u_i\)與 \(\alpha = \frac{p(x)}{k\cdot q(x)}\):
如果 \(u_i\le \alpha\) 則認為 \(x_i\)為服從p(x)的有效樣本，反之，則認為無效，舍棄。
重復上述步驟直至抽樣結束。其中\(\alpha\)稱為接受率，這是此方法名字的由來。
下圖為實驗結果：
接受拒絕采樣的證明，即是證明，經過上述過程得到的樣本是否服從目標分布p(x)：
設目標分布與建議分布分別用T，Q表示,其對應的概率密度函數分別為t(x)與q(x), 證明此采樣方法的有效性即是證明：
其中，\(U\)為服從(0,1)均勻分布的隨機變量。
蒙特卡洛方法
蒙特卡洛方是一個賭場的名字，其作為隨機模擬方法的名稱有說是因為玩笑也有人說是因為機密。
蒙特卡洛方法最初目的是解決直接求解積分（或求和）有困難的問題，因為積分不能直接求得，需要通過某種方法來近似，蒙特卡洛方法即采用隨機抽樣的方法，在大數定律的作用下，隨着抽樣次數增加，近似會越來越精確。
設在[a,b]上，函數\(0 \le f(x)\le M\), 計算：
圖中陰影部分被矩形所包，欲求此部分面積，可在這個矩形內隨機的投點，某點落在陰影部分的概率為:
根據大數定律，當投點次數很大時，可將 \(\frac{n}{N}\)作為P的近似值，其中n為落入陰影中的點數，而N為投點總數。於是：
如果知道\(f(x)\)，可以在[a,b]上隨機的抽樣\(\{x_0,x_2,...,x_{n-1}\}\). 然后計算分別計算\(，求得陰影面積f(x_i)，求得陰影面積\) :
重要性采樣：
重要性采樣方法作為一種蒙特卡洛方法，可降低估計方差。
重要性采樣考慮的是通過蒙特卡洛方法求取\(\mu = E[f(X)]\)， 其中X是服從分布\(p(x)\)：
當\(p(x)\)不容易抽樣時，重要性采樣如接受-拒絕采樣一樣引入了另一個相對容易采樣的分布，不過重要性采樣沒有對樣本進行接受與否的判定，而是對所有樣本進行加權平均，其想法也是非常的直觀:
其中， \(w(x) = \frac{p(x)}{q(x)}\).稱為重要性采樣比率（Importance Sampling Ratio 或 likelihood ratio）.
MCMC(Markov Chain Monte Carlo)
在隨機過程中，馬爾可夫性質表示當給定t時刻的狀態，那么t+1時刻隨機過程的狀態概率就是確定的：
即狀態轉移的概率只依賴於前一個狀態，而與更早的狀態無關。
舉個例子，某生物實驗室，為科研需要培養大量的斑馬魚，未經實驗的斑馬魚假設只有紅黃藍三種顏色，實驗室的科研人員發現在穩定的實驗室環境內，三種顏色的魚都會生出其色顏色的斑馬魚（為簡單，這里假定只有相同顏色的斑馬魚才會交配。），比如紅色的斑馬魚其后代有0.72概率還是紅色，而有0.17的概率生成黃色斑馬魚，有0.11的概率生成藍色斑馬魚。具體的：
發現這了這個轉移概率矩陣后，科研人員對每一代的斑馬魚的比例都了若指掌。比如有一次他們進購了三種魚的比率是【0.2，0.3，0.5】.
| Generation | Red | Blue | Yellow |
|---|---|---|---|
| 0 | 0.2 | 0.3 | 0.5 |
| 1 | 0.303 | 0.461 | 0.236 |
| 2 | 0.399 | 0.298 | 0.304 |
| 3 | 0.422 | 0.343 | 0.235 |
| 4 | 0.445 | 0.301 | 0.254 |
| 5 | 0.45 | 0.313 | 0.236 |
| 6 | 0.456 | 0.302 | 0.242 |
| 7 | 0.457 | 0.306 | 0.237 |
| 8 | 0.458 | 0.303 | 0.238 |
| 9 | 0.459 | 0.304 | 0.237 |
| 10 | 0.459 | 0.303 | 0.238 |
| 11 | 0.459 | 0.304 | 0.237 |
| 12 | 0.459 | 0.303 | 0.237 |
| 13 | 0.459 | 0.303 | 0.237 |
| 14 | 0.459 | 0.303 | 0.237 |
| ... | ... | ... | ... |
當到第11代時，三種顏色的比例就穩定了。科研人員覺得很有意思，於是他們又一次購買魚苗時，將比例變成了[0.05,0.94,0.01] 這樣一個比較奇怪的比例：
| Generation | Red | Blue | Yellow |
|---|---|---|---|
| 0 | 0.05 | 0.94 | 0.01 |
| 1 | 0.347 | 0.148 | 0.505 |
| 2 | 0.359 | 0.468 | 0.172 |
| 3 | 0.434 | 0.259 | 0.307 |
| 4 | 0.435 | 0.346 | 0.219 |
| 5 | 0.454 | 0.291 | 0.255 |
| 6 | 0.453 | 0.315 | 0.232 |
| 7 | 0.458 | 0.3 | 0.242 |
| 8 | 0.458 | 0.306 | 0.236 |
| 9 | 0.459 | 0.302 | 0.239 |
| 10 | 0.459 | 0.304 | 0.237 |
| 11 | 0.459 | 0.303 | 0.238 |
| 12 | 0.459 | 0.304 | 0.237 |
| 13 | 0.459 | 0.303 | 0.237 |
| 14 | 0.459 | 0.303 | 0.237 |
| ... | ... | ... | ... |
同樣地，經過幾代之后，比例再次穩定且與上次相同。這不是巧合，這是絕大多數馬氏鏈都具有的收斂性質：
對於一個非周期的馬爾可夫鏈，其狀態轉移概率矩陣為P，且其各個狀態都是連通的，那么\(\lim _{n \rightarrow \infty} P_{i j}^{n}\)存在且與 i無關，如果記\(\lim _{n \rightarrow \infty} P_{i j}^{n}=\pi(j)\),那么我們可得到：
2）
- \(\pi\)是 方程\(\pi P=\pi\)唯一非負解，其中\(\pi=[\pi(1), \pi(2), \cdots, \pi(j), \cdots], \quad \sum_{i=0}^{\infty} \pi_{i}=1\)，稱為馬氏鏈的平穩分布。（PS: \(\pi\) 是概率分布，不是某一確定數值。）
這個是馬爾可夫鏈很重要的收斂性質，后面的所有的都是基於此性質的。在這個性質中，'非周期'表示的是隨機過程的狀態轉化不是循環的，即對任意狀態i，d為集合\(\left\{n | n \geq 1, P_{i i}^{n}>0\right\}\)的最大公約數，其值為1則表明馬氏鏈是非周期的。 '各個狀態連通'表明狀態轉移概率矩陣中各個元素都是非0的。對2)中有證明也非常的直觀：
同時對兩邊取極限，即可得到\(\pi(j)=\sum_{i=0}^{\infty} \pi(i) P_{i j}\). 另外，馬爾可夫鏈的狀態數可是有限的，也可以是無限的，即上述性質適用於離散概率分布與連續概率分布。這里所謂“平穩分布”即是指經概率轉移矩陣的作用，分布不再變化。
當已知某平穩分布對應的轉移概率矩陣P，就可得到此平穩分布的樣本集。由上述性質可知，無論初始概率分布（\(\pi_0(x)\)）為何，經過多輪（比如n輪）概率轉移后，就會收斂到平穩概率分布，即:
因此，可從簡單分布中抽樣，得到初始樣三本\(x_0\)，然后根據條件概率分布\(P(x_1|x_0)\)得到\(x_1\), 依次類推，根據\(P(x_2|x_1),P(x_3|x_2),...,P(x_n|x_{n-1}),...,P(x_{n+m}|p_{n+m-1})\) 分別得到\(x_2,x_3,...,x_n,...,x_{n+m}\),其中\((x_{n+1},x_{n+2},...,x_{n_m})\) 即是服從所需平穩分布的樣本集。
上面的收斂性質給予我們很大的啟發，對於某分布\(\pi(x)\)，要得到其對應的樣本，只需構造一個具有轉移概率矩陣為P的馬爾可夫鏈，使其收斂的平穩分布為\(p(x)\)。那么如何做呢？這就涉及一個很重要的概念，稱為細致平衡條件：如果非周期馬爾可夫鏈的狀態轉移矩陣P的概率分布\(\pi(x)\)對於所有的i,j滿足：
則稱為概率分布\(\pi(x)\)是狀態轉移矩陣P的平穩分布。
這個條件的很直觀：
不過，通常情況下，這個條件無法滿足，即很難找到這樣一個馬爾可夫鏈對應的轉移概率矩陣，使得上式成立。比如轉移概率為Q的馬爾可夫鏈，可以驗證，一般地：
為使等式成立，可以引入另一個矩陣\(\alpha\)，即：
為達此目的，\(\alpha\)需要滿足些條件：
令：
\(\alpha\)稱為接受率,類似上面提到的接受-拒絕采樣中的接受率，在接受-拒絕采樣中，接受率調節的是從易抽樣分布中得到的樣本是否接受其成為目標分布的樣本，而在這里，接受率調節服從Q的狀態轉移是否接受成為P的狀態轉移。這里的狀態轉移矩陣P就是分布\(\pi(x)\)對應的馬爾可夫鏈的轉移概率矩陣P。即滿足細致平穩條的馬爾可夫鏈的轉移概率矩陣P可通過接受率的調節，從一個具有普通的轉移概率矩陣Q的馬爾可夫鏈變換得到。
這樣，就得到了MCMC采樣的一般過程：
1） 確定目標分布\(\pi(x)\),
2）設定收斂前的轉移次數\(n_1\), 所需樣本數\(n_2\). 選定一個馬爾可夫鏈的概率轉移矩陣Q，只需馬爾可夫鏈非周期，各狀態之間相互連通即可。
3）從簡單分布如均勻分布中采樣得到初始樣本\(x_0\),
- loop: i from 0 到 \(n_1+n_2 -1\):
- 根據條件概率\(Q(x|x_i)\) 得到 x';
- 從[0,1]的均勻分布中抽樣得到 \(\mu\);
- 如果 \(\mu \lt \alpha(x_i,x') = \pi(x')Q(x',x_i)\)：
- \(x_{i+1} = x'\),即接受 \(x_i\)到\(x'\)的轉移，
- 否則,不接受轉移：
- \(x_{i+1} = x_i\)
- 收集樣本\(\{x_{n_1},x_{n_1+1},...,.x_{n_1+n_2}\}\),即為平稱分布\(\pi(x)\)的樣本集。
此方法對離散、連續情形均適用。
可見，上述方法非常的強大，適用范圍很廣，但實際應用過程中發現，上述過程收斂很慢，要將\(n_1\)設置的非常大，且設定比較困難，而且收斂過慢，不但耗時且浪費計算資源。這個問題由Metropolis 與 Hastings 解決，並產生了Metropolis-Hastings算法（或稱M-H算法）。
M-H 算法
MCMC采樣效率低的原因是接受率通常比較小，導致大部分狀態轉移不被接受。比如上面過程中\(\alpha(x_i,x' ) = 0.1\),那么\(\mu\lt \alpha(x_i,x')\)的概率就只有0.1。優化上述問題的關鍵在於提高接受率。在細致平穩條件中：
兩邊同時乘一個系數如 k，等式仍是成立的：
如果令:
成為新的接受率，且\(k\gt 1\), 那么這個接受率在不違背細致平穩條件的情況下被放大了！因為接受率不大於1，放大接受率的上限為1，那不妨讓\(\alpha'(i,j), \alpha'(j,i)\) 中較大的為1，這樣算法抽樣的效率達到最高。即：
這樣，就得到了M-H采樣算法的一般過程,與MCMC基本一致，只是\(\mu\)的接受空間變大了：
1） 確定目標分布\(\pi(x)\),
2）設定收斂前的轉移次數\(n_1\), 所需樣本數\(n_2\). 選定一個馬爾可夫鏈的概率轉移矩陣Q，只需馬爾可夫鏈非周期，各狀態之間相互連通即可。
3）從簡單分布如均勻分布中采樣得到初始樣本\(x_0\),
- loop: i from 0 到 \(n_1+n_2 -1\):
- 根據條件概率\(Q(x|x_i)\) 得到 x';
- 從[0,1]的均勻分布中抽樣得到 \(\mu\);
- 如果 \(\mu \lt \alpha(x_i,x') =\min \left\{\frac{\pi(x') Q(x', x_i)}{\pi(x_i) Q(x_i, x')}, 1\right\}\)：
- \(x_{i+1} = x'\),即接受 \(x_i\)到\(x'\)的轉移，
- 否則,不接受轉移：
- \(x_{i+1} = x_i\)
- 收集樣本\(\{x_{n_1},x_{n_1+1},...,.x_{n_1+n_2}\}\),即為平稱分布\(\pi(x)\)的樣本集。
Gibbs Sampling
對於高維數據，上述的M-H算法再次突顯出了弊端。對於高維數據，聯兩數據的聯合分布會變得很復雜，甚至無法求解；另外，接受率的存在，計算接受率，並在其限制下接受樣本，都使得很多算力浪費，效率不高。於是Gibbs進一步提出了Gibbs采樣算法。當維度大於1時，不妨先從二維數據開始探討。設兩個點\(A(x_1,y_1),B(x_1,y_2)\), 服從概率分布\(p(x)\). 這兩個點的第一維度坐標相同，可推出：
進而：
也就是：
可見當其x維度相同時，兩個點在條件概率\(p(y|x)\)作為轉移概率的情況下是滿足細致平穩條件的。同理，如果有一點\(C(x_2,y_1)\), 那么點\(A,C\)滿足細致平穩條件：
即：
因此對於二維數據的Gibbs采樣：
-
確定目標分布 \(\pi(x_1,x_2)\), 並狀態轉移次數\(n_1\)與樣本個數\(n_2\);
-
隨機初始化狀態值\(x_1^0,x_2^0\),
-
loop(0， \(n_1+n_2\)):
- 從條件概率分布\(P(x_2|x_1^t)\)中采樣得到\(x_2^{t+1}\);
- 從條件概率分布P(x_1|x_2^{t+1})中采樣得到樣本\(x_1^{t+1}\)
-
上述循環得到的\(\left\{\left(x_{1}^{\left(n_{1}\right)}, x_{2}^{\left(n_{1}\right)}\right),\left(x_{1}^{\left(n_{1}+1\right)}, x_{2}^{\left(n_{1}+1\right)}\right), \dots,\left(x_{1}^{\left(n_{1}+n_{2}-1\right)}, x_{2}^{\left(n_{1}+n_{2}-1\right)}\right)\right\}\) 樣本即為所需樣本。
這個公式很容易推廣到 n 維數據，在n-1個維度固定的情況下，所選維度上的轉換是服從細致平穩條件的，n維Gibbs采樣過程：
-
確定目標分布 \(\pi(x_1,x_2, ...,x_n)\), 並狀態轉移次數\(n_1\)與樣本個數\(n_2\);
-
隨機初始化狀態值\(x_1^0,x_2^0,...,x_n^0\)
-
loop(0， \(n_1+n_2\)):
- 從條件概率分布\(P(x_1|x_2^t,x_3^t,...,x_n^t)\)中采樣得到$x_1^{t+1} $;
- 從條件概率分布\(P(x_2|x_1^{t+1},x_3^t,...,x_n^t)\)中采樣得到樣本$x_2^{t+1} $;
- ...
- 從條件概率分布\(P(x_j|x_1^{t+1},x_2^{t+1},...,x_{j-1}^{t+1},x_{j+1}^t,...,x_n^t)\)中采樣得到樣本\(x_j^{t+1}\);
- ...
- 從條件概率分布\(P(x_n|x_1^{t+1},x_2^{t+1},...,x_{n-1}^{t+1})\)中采樣得到樣本\(x_n^{t+1}\).
-
上述循環得到的\(\left\{\left(x_{1}^{\left(n_{1}\right)}, x_{2}^{\left(n_{1}\right)}, \dots, x_{n}^{\left(n_{1}\right)}\right), \ldots,\left(x_{1}^{\left(n_{1}+n_{2}-1\right)}, x_{2}^{\left(n_{1}+n_{2}-1\right)}, \ldots, x_{n}^{\left(n_{1}+n_{2}-1\right)}\right)\right\}\) 樣本即為所需樣本。