Hulu機器學習問題與解答系列 | 十四：如何對高斯分布進行採樣

「如何對高斯分布進行採樣」

[場景描述]

高斯分布，又稱正態分布，是一個在數學、物理及工程領域都非常重要的概率分布。在實際應用中，我們經常需要對高斯分布進行採樣。雖然在很多編程語言中，直接調用一個函數就可以生成高斯分布隨機數，但了解其中的具體演算法能夠加深我們對相關概率統計知識的理解；此外，高斯分布的採樣方法有多種，通過展示不同的採樣方法在高斯分布上的具體操作以及性能對比，我們會對這些採樣方法有更直觀的印象。

[問題描述]

如果讓你來實現一個高斯分布隨機數生成器，你會怎麼做？

背景知識：概率統計

[解答與分析]

首先，假設隨機變數z服從標準正態分布N(0,1)，令

x = σ·z + μ

則x服從均值為μ、方差為σ2的高斯分布N(μ, σ2)。因此，任意高斯分布都可以由標準正態分布通過拉伸和平移得到，所以這裡我們只考慮標準正態分布的採樣。另外，幾乎所有的採樣方法都是以均勻分布隨機數作為基本操作，因此這裡假設我們已經有均勻分布隨機數生成器了。均勻分布隨機數一般用線性同餘法來生成（偽隨機數），具體參見文獻[1]。

常見的採樣方法有逆變換法 (Inverse Transform Method)、拒絕採樣法 (Rejection Sampling)、重要性採樣及其重採樣 (Importance Sampling, Sampling-Importance-Resampling)、馬爾科夫蒙特卡洛採樣法 (Markov Chain Monte Carlo) 等。具體到高斯分布，我們需要如何採樣呢？

如果直接用逆變換法，基本操作如下：

上述逆變換法需要求解erf(x)的逆函數，這並不是一個初等函數，沒有顯式解，計算起來比較麻煩。為了避免這種非初等函數的求逆操作，Box-Muller演算法採用如下解決方案：既然單個高斯分布的累計分布函數不好求逆，那麼兩個獨立的高斯分布的聯合分布呢？假設x, y是兩個服從標準正態分布的獨立隨機變數，它們的聯合概率密度為：

考慮(x, y)在圓盤{(x, y) | x2 + y2 ≤ R2}上的概率：

通過極坐標變換將(x, y)轉化為(r, θ)，可以很容易求得上述二重積分：

這裡F(R)可以看成是極坐標中r的累積分布函數。由於F(R)的計算公式比較簡單，逆函數也很容易求得，所以可以利用逆變換法來對r進行採樣；對於θ，在[0, 2π]上進行均勻採樣即可。這樣就得到了(r,θ)，經過坐標變化即可得到符合標準正態分布的(x, y)。具體採樣過程如下：

Box–Muller演算法由於需要計算三角函數，相對來說還是比較耗時，而Marsaglia polar method則避開了三角函數的計算，因而更快，其具體採樣操作如下：

除了逆變化法，我們還可以利用拒絕採樣法，選擇一個比較好計算累積分布逆函數的參考分布來覆蓋當前正態分布（可以乘以一個常數倍），進而轉化為對參考分布的採樣以及對樣本點的拒絕/接收操作。考慮到高斯分布的特性，這裡可以用指數分布來作為參考分布。指數分布的累積分布及其逆函數都比較容易求解。由於指數分布的樣本空間為x≥0，而標準正態分布的樣本空間為(-∞, +∞)，因此還需要利用正態分布的對稱性來在半坐標軸和全坐標軸之間轉化。具體來說，取λ=1的指數分布作為參考分布，

實際應用時，M需要儘可能小，這樣每次的接受概率大，採樣效率更高。因此，可以取

因此，具體的採樣過程如下：

拒絕採樣法的效率取決於接受概率的大小：參考分布與目標分布越接近，則採樣效率越高。有沒有更高效的拒絕採樣演算法呢？這就是Ziggurat演算法，該演算法本質也是拒絕採樣，但採用多個階梯矩形來逼近目標分布（如圖1所示）。Ziggurat演算法雖然看起來稍微繁瑣，但實現起來並不複雜，操作也非常高效，其具體採樣步驟可以參考文獻[2]。

[總結與擴展]

高斯分布的採樣方法還有很多，我們只列舉了幾種最常見的。具體面試時，候選人不需要回答所有的方法，知道其中一兩種即可，面試官可以針對這一兩種方法深入提問，如理論證明、優缺點、性能等。如果候選人沒有思路，面試官可以引導其回憶那些通用的採樣方法，如何將那些策略用到高斯分布這個具體案例上。另外，本題還可以適當擴展，把一般的高斯分布換成截尾高斯分布 (Truncated Gaussian Distribution) ，如何採樣？如果是高維隨機變數，拒絕採樣法會存在什麼問題？怎麼解決呢？

[參考文獻]

[1] Linear congruential generator,

Linear congruential generator?

en.wikipedia.org

[2] Ziggurat algorithm,