概率與數理統計-蒙特卡洛之拒絕采樣

聲明

轉：> https://blog.csdn.net/u010159842/article/details/78959515

介紹

蒙特卡洛（Monte Carlo）方法是二十世紀四十年代中期由於科學技術的發展和電子計算機的發明，而被提出的一種以概率統計理論為基礎的數值計算方法。它的核心思想就是使用隨機數（或更常見的偽隨機數）來解決一些復雜的計算問題。

當所求解問題可以轉化為某種隨機分布的特征數（比如隨機事件出現的概率，或者隨機變量的期望值等）時，往往就可以考慮使用蒙特卡洛方法。通過隨機抽樣的方法，以隨機事件出現的頻率估計其概率，或者以抽樣的數字特征估算隨機變量的數字特征，並將其作為問題的解。這種方法多用於求解復雜的高維積分問題。

實際應用中，我們所要面對的第一個問題就是如何抽樣？在統計學中， 抽樣（或稱采樣）是指從目標總體中抽取一部分個體作為樣本的過程。

例如，我們想知道一所大學里所有男生的平均身高。但是因為學校里的男生可能有上萬人之多，所以為每個人都測量一下身高可能存在困難，於是我們從每個學院隨機挑選出100名男生來作為樣本，這個過程就是抽樣

但是在計算機模擬時，我們所說的抽樣，其實是指從一個概率分布中生成觀察值（observations）的方法。而這個分布通常是由其概率密度函數（PDF）來表示的。而且，即使在已知PDF的情況下，讓計算機自動生成觀測值也不是一件容易的事情。從本質上來說，計算機只能實現對均勻分布（Uniform distribution）的采樣。

具體來說，我們可能要面對的問題包括：

• 計算機只能實現對均勻分布的采樣，但我們仍然可以在此基礎上對更為復雜的分布進行采樣，那具體該如何操作呢？
• 隨機分布的某些數字特征可能需要通過積分的形式來求解，但是某些積分可能沒有（或者很難求得）解析解，彼時我們該如何處理呢？
• 在貝葉斯推斷中，后驗概率的分布是正⽐於先驗和似然函數之積的，但是先驗和似然函數的乘積形式可能相對復雜，我們又該如何對這種形式復雜的分布進行采樣呢？

Inverse CDF 方法

比較簡單的一種情況是，我們可以通過PDF與CDF之間的關系，求出相應的CDF。或者我們根本就不知道PDF，但是知道CDF。此時就可以使用Inverse CDF的方法來進行采樣。這種方法又稱為逆變換采樣（Inverse transform sampling）。

如果你對PDF和CDF的概念有點模糊，我們不妨先來一起回顧一下它們的定義。對於隨機變量\(X\)，如下定義的函數 \(F\):

\[F(x)=P\{X<= x\}, -\infty < x < \infty \]

稱為 \(X\) 的累積分布函數（CDF，Cumulative Distribution Function）。對於連續型隨機變量 \(X\) 的累積分布函數 \(F(x)\)，如果存在一個定義在實數軸上的非負函數 \(f(x)\)，使得對於任意實數 \(x\)，有下式成立：

\[F(x)=\int_{-\infty}^{\infty}{f(t)}dt \]

則稱 \(f(x)\) 為 \(X\) 的概率密度函數（PDF，Probability Density Function）。顯然，當概率密度函數存在的時候，累積分布函數是概率密度函數的積分。

所以，通常我們可以通過對PDF（如下圖中的左圖所示為正態分布的PDF）進行積分來得到概率分布的CDF（如下圖中的右圖所示為正態分布的CDF）。然后我們再得到CDF的反函數 \(F^{-1}(u)\)，如果你想得到 \(m\) 個觀察值，則重復下面的步驟 \(m\) 次:

• 從 Uniform(0,1) 中隨機生成一個值（前面已經說過，計算機可以實現從均勻分布中采樣），用 \(u\) 表示。
• 計算 \(F_{-1}(u)\) 的值 \(x\)，則 \(x\) 就是從 \(f(x)\) 中得出的一個采樣點。

以下圖為例，如果從 Uniform(0,1) 中隨機生成的值 \(u=0.8413\)，則可以算得\(F_{-1}(u)=1\)，則此次從正態分布中生成的隨機數就是 1。

你可能會好奇，面對一個具有復雜表達式的函數， Inverse CDF 方法真的有效嗎？來看下面這個例子。假設現在我們希望從下面這個PDF中抽樣：

\[f(x)=\left\{\begin{array}{ll}{8 x} & {, \text {if } 0 \leq x<0.25} \\ {\frac{8}{3}-\frac{8}{3} x} & {, \text { if } 0.25 \leq x \leq 1} \\ {0} & {, \text { otherwise }}\end{array}\right. \]

可以算得相應的CDF為

\[F(x)=\left\{\begin{array}{ll}{0} & {, \text { if } x<0} \\ {4 x^{2}} & {, \text { if } 0 \leq x<0.25} \\ {\frac{8}{3} x-\frac{4}{3} x^{2}-\frac{1}{3}} & {, \text { if } 0.25 \leq x \leq 1} \\ {1} & {, \text { if } x>1}\end{array}\right. \]

對於 \(u∈[0,1]\)，它的反函數為:

\[F^{-1}(u)=\left\{\begin{array}{ll}{\frac{\sqrt{u}}{2}} & {, \text { if } 0 \leq u<0.25} \\ {1-\frac{\sqrt{3(1-u)}}{2}} & {, i f 0.25 \leq u \leq 1}\end{array}\right. \]

拒絕采樣

我們已經看到 Inverse CDF 方法確實有效。但其實它的缺點也是很明顯的，那就是有些分布的 CDF 可能很難通過對 PDF 的積分得到，再或者 CDF 的反函數也很不容易求。這時我們可能需要用到另外一種采樣方法，這就是我們即將要介紹的拒絕采樣。

下面這張圖很好地闡釋了拒絕采樣的基本思想。假設我們想對 PDF 為 \(p(x)\) 的函數進行采樣，但是由於種種原因（例如這個函數很復雜），對其進行采樣是相對困難的。但是另外一個 PDF 為 \(q(x)\) 的函數則相對容易采樣，例如采用 Inverse CDF 方法可以很容易對對它進行采樣，甚至 \(q(x)\) 就是一個均勻分布（別忘了計算機可以直接進行采樣的分布就只有均勻分布）。那么，當我們將 \(q(x)\) 與一個常數 MM 相乘之后，可以實現下圖所示之關系，即\(M⋅q(x)\) 將 \(p(x)\) 完全“罩住”

你當然可以采用嚴密的數學推導來證明Reject Sampling的可行性。但它的原理從直觀上來解釋也是相當容易理解的。你可以想象一下在上圖的例子中，從哪些位置抽出的點會比較容易被接受。顯然，紅色曲線和綠色曲線所示之函數更加接近的地方接受概率較高，也即是更容易被接受，所以在這樣的地方采到的點就會比較多，而在接受概率較低（即兩個函數差距較大）的地方采到的點就會比較少，這也就保證了這個方法的有效性。

Reject Sampling方法確實可以解決我們的問題。但是它的一個不足涉及到其采樣效率的問題。
最理想的情況下，參考分布應該跟目標分布越接近越好，從圖形上來看就是包裹的越緊實越好。但是這種情況的參考分布往往又不那么容易得到。在滿足某些條件的時候也確實可以采用所謂的改進方法，即Adaptive Rejection Sampling。

自適應的拒絕采樣（Adaptive Rejection Sampling）

前面我們已經分析了，拒絕采樣的弱點在於當被拒絕的點很多時，采樣的效率會非常不理想。同時我們也支持，如果能夠找到一個跟目標分布函數非常接近的參考函數，那么就可以保證被接受的點占大多數（被拒絕的點很少）。這樣一來便克服了拒絕采樣效率不高的弱點。如果函數是 log-concave 的話，那么我們就可以采樣自適應的拒絕采樣方法。什么是 log-concave 呢？還是回到我們之前介紹過的 Beta 分布的PDF，我們用下面的代碼來繪制 Beta(2, 3) 的函數圖像，以及將 Beta(2, 3) 的函數取對數之后的圖形

(......更詳細的內容看原文)