3.2　蒙特卡洛采樣

在3.1節中，通過求解定積分這個具體的例子來演示了蒙特卡洛方法的基本思想，而其中的核心就是使用隨機數。當所求解問題可以轉化為某種隨機分布的特征數（如隨機事件出現的概率，或者隨機變量的期望值等）時，往往就可以考慮使用蒙特卡洛方法。通過隨機抽樣的方法，以隨機事件出現的頻率估計其概率，或者以抽樣的數字特征估算隨機變量的數字特征，并將其作為問題的解。這種方法多用于求解復雜的高維積分問題。

實際應用中，所要面對的第一個問題就是如何抽樣？注意，在計算機模擬時，這里所說的抽樣其實是指從一個概率分布中生成觀察值（observations）的方法。而這個分布通常是由其概率密度函數來表示的。前面曾經提過，即使在已知PDF的情況下，讓計算機自動生成觀測值也不是一件容易的事情。從本質上來說，計算機只能實現對均勻分布的采樣。幸運的是，仍然可以在此基礎上對更為復雜的分布進行采樣。

3.2.1　逆采樣

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

假設已經得到了CDF的反函數F^-1（u），如果想得到m個觀察值，則重復下面的步驟m次：

（1）從U（0，1）中隨機生成一個值（計算機可以實現從均勻分布中采樣），用u表示；

（2）計算F^-1（u）的值x，則x就是從目標分布f（x）中得出的一個采樣點。

面對一個具有復雜表達式的函數，逆變換采樣法真有效嗎？來看一個例子，假設現在希望從具有下面這個PDF的分布中采樣

可以算得相應的CDF為

對于u∈［0，1］，上述CDF的反函數為

下面在R中利用上面的方法采樣10000個點，并以此來演示抽樣的效果。

將所得結果與真實的PDF函數圖形進行對照，如圖3-5所示。可見由逆變換采樣法得到的點所呈現處理的分布與目標分布非常吻合。

下面再舉一個稍微復雜一點的例子，已知分布的PDF如下

可以算得相應的CDF為

對于u∈［0，1］，它的反函數為

同樣，給出R中的示例代碼如下。

下面這段代碼利用R中提供的一些內置函數實現了已知PDF時基于逆變換方法的采樣，將新定義的函數命名為samplepdf（）。當然，對于那些過于復雜的PDF函數（例如很難積分的），samplepdf（）確實有力所不及的情況。但是對于標準的常規PDF，該函數的效果還是不錯的。

下面就用samplepdf（）函數來對上面給定的h（x）進行采樣，然后再跟之前所得結果進行對比。

代碼執行結果如圖3-6所示。

3.2.2　博克斯-穆勒變換

博克斯-穆勒變換（Box-Muller Transform）最初由喬治·博克斯（George Box）與默文·穆勒（Mervin Muller）在1958年共同提出。博克斯是統計學的一代大師，統計學中的很多名詞術語都以其名字命名。博克斯與統計學的家學淵源相當深厚，他的導師是統計學開山鼻祖皮爾遜的兒子，英國統計學家埃貢·皮爾遜（Egon Pearson），博克斯還是統計學的另外一位巨擘級奠基人費希爾的女婿。統計學中的名言“所有模型都是錯的，但其中一些是有用的”也出自博克斯之口。

本質上來說，計算機只能生產符合均勻分布的采樣。如果要生成其他分布的采樣，就需要借助一些技巧性的方法。而在眾多的“其他分布”中，正態分布無疑占據著相當重要的地位。下面這個定理，就為生成符合正態分布的采樣（隨機數）提供了一種方法，而且這也是很多軟件或者編程語言的庫函數中生成正態分布隨機數時所采樣的方法。

定理（Box-Muller變換）：如果隨機變量U₁和U₂是獨立同分布的，且U₁，U₂~U［0，1］，則

其中，Z₀和Z₁獨立且服從標準正態分布。

如何來證明這個定理呢？這需要用到一些微積分中的知識，首先回憶一下二重積分化為極坐標下累次積分的方法

假設現在有兩個獨立的標準正態分布X~N（0，1）和Y~N（0，1），由于兩者相互獨立，則聯合概率密度函數為

做極坐標變換，則x=Rcosθ，y=Rsinθ，則有

這個結果可以看成是兩個概率分布的密度函數的乘積，其中一個可以看成是［0，2π］上均勻分布，將其轉換為標準均勻分布則有θ~U（0，2π）=2πU₂。

另外一個的密度函數為

則其累計分布函數CDF為

這個CDF函數的反函數可以寫成

根據逆變換采樣的原理，如果有個PDF為P（R）的分布，那么對齊CDF的反函數進行均勻采樣所得的樣本分布將符合P（R）的分布，而如果u是均勻分布的，那么U₁=1-u也將是均勻分布的，于是用U₁替換1-u，最后可得

結論得證。最后來總結一下利用Box-Muller變換生成符合高斯分布的隨機數的方法：

（1）產生兩個隨機數U₁，U₂~U［0，1］；

（2）用它們來創造半徑和和夾角θ=2πU₂；

（3）將（R，θ）從極坐標轉換到笛卡兒坐標：（Rcosθ，Rsinθ）。

3.2.3　拒絕采樣與自適應拒絕采樣

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

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

然后重復如下步驟，直到獲得m個被接受的采樣點：

（1）從q（x）中獲得一個隨機采樣點x_i；

（2）對于x_i計算接受概率（acceptanceprobability）

（3）從U（0，1）中隨機生成一個值，用u表示；

（4）如果α≥u，則接受x_i作為一個來自p（x）的采樣值，否則就拒絕x_i并回到第一步。

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

還是以本章前面給出的那個分段函數f（x）為例來演示拒絕采樣方法。如圖3-8所示，所選擇的參考分布是均勻分布（當然也可以選擇其他的分布，但采用均勻分布顯然是此處最簡單的一種處理方式）。而且令常數M=3。

下面給出R中的示例代碼，易見此處采樣點數目為10 000。

上述代碼的執行結果如圖3-9所示，可見采樣結果是非常理想的。

下面的例子演示了對（表達式非常復雜的）beta（3，6）分布進行拒絕采樣的效果。這里采用均勻分布作為參考分布。而且這里的Mq（x）所取之值就是beta（3，6）分布的極大值，它的函數圖形應該是與beta（3，6）的極值點相切的一條水平直線。

圖3-10給出了采樣50 000個點后的密度分布情況，可見采樣分布與目標分布beta（3，6）非常吻合。

拒絕采樣的方法確實可以解決我們的問題。但是它的一個不足涉及其采樣效率的問題。針對上面給出的例子而言，我們選擇了離目標函數最近的參考函數，就均勻分布而言，已經不能有更進一步的方法了。但即使這種，在這個類似鐘形的圖形兩側其實仍然會拒絕掉很多很多采樣點，這種開銷相當浪費。最理想的情況下，參考分布應該跟目標分布越接近越好，從圖形上來看就是包裹的越緊實越好。但是這種情況的參考分布往往又不那么容易得到。在滿足某些條件的時候也確實可以采用所謂的改進方法，即自適應的拒絕采樣（Adaptive Rejection Sampling）。

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

上述代碼的執行結果如圖3-11所示，其中（a）是beta（2，3）的概率密度函數圖形，（b）是將beta（2，3）的函數取對數之后的圖形，你可以發現結果是一個凹函數（concave）。那么beta（2，3）就滿足log-concave的要求。

然后在對數圖像上找一些點做圖像的切線，如圖3-12所示。因為取對數后的函數是凹函數，所以每個切線都相當于一個超平面，而且對數圖像只會位于超平面的一側。

同時給出用以繪制圖3-12的代碼，要知道R語言的一個強項就是繪圖。

再把這些切線轉換回原始的beta（2，3）圖像中，顯然原來的線性函數會變成指數函數，它們將對應圖3-13中的一些曲線，這些曲線會被原函數的圖形緊緊包裹住。特別是當這些的指數函數變得很多很稠密時，以彼此的交點作為分界線，其實相當于得到了一個分段函數。這個分段函數是原函數的一個逼近。用這個分段函數來作為參考函數再執行拒絕采樣，自然就完美地解決了之前的問題。

下面是用來畫出圖3-13的R語言代碼。

這無疑是一種絕妙的想法。而且這種想法，在前面其實已經暗示過。在上一部分最后一個例子中，我們其實就是選擇了一個與原函數相切的均勻分布函數來作為參考函數。我們當然會想去選擇更多與原函數相切的函數，然后用這個函數的集合來作為新的參考函數。只是由于原函數的凹凸性無法保證，所以直線并不是一種好的選擇。而自適應拒絕采樣（Adaptive Rejection Sampling，ARS）所采用的策略則非常巧妙地解決了我們的問題。當然函數是log-concave的條件必須滿足，否則就不能使用ARS。

下面給出一個在R中進行自適應拒絕采樣的例子。顯然，該例子要比之前的代碼簡單許多。因為R中ars包已經提供了一個現成的用于執行自適應拒絕采樣的函數，即ars（）。關于這個函數在用法上的一些細節，讀者還可以進一步參閱R的幫助文檔，這里不再贅言。此次我們需要指出：ars（）函數中兩個重要參數，一個是對原分布的PDF取對數，另外一個則是對PDF的對數形式再進行求導（在求導時我們忽略了前面的系數項），其實也就是為了確定切線。

上述代碼的執行結果如圖3-14所示。