詳解Box-Muller方法生成正態(tài)分布

本公眾號MyEncyclopedia定期發(fā)布AI，算法，工程，大數(shù)據(jù)交叉領(lǐng)域的深度和前沿文章。歡迎關(guān)注，收藏和點贊。公眾號內(nèi)有本文對應(yīng)的配套的視頻講解。

在學(xué)習(xí)了一些基本的統(tǒng)計變量生成法之后，這次我們來看看如何生成正態(tài)分布。它就是大名鼎鼎的 Box-Muller 方法，Box-Muller 的理解過程可以體會到統(tǒng)計模擬的一些精妙思想。

嘗試逆變換方法

關(guān)于逆變換方法，在用逆變換采樣方法構(gòu)建隨機變量生成器中有詳細的講解，那么我們就先嘗試通過逆變換方法標(biāo)準流程來生成正態(tài)分布。

正態(tài)分布的 PDF 表達式為

對應(yīng)的函數(shù)圖形是鐘形曲線

根據(jù) PDF，其 CDF 的積分形式為

和所有 PDF CDF 關(guān)系一樣， 表示 累積到 點的面積。

很不幸的是， 無法寫出一般數(shù)學(xué)表達式，因而也無法直接用逆變換方法。

二維映射到一維

我們知道，高維正態(tài)分布有特殊的性質(zhì)：它的每一維的分量都是正態(tài)分布；單個維度對于其他維度的條件概率分布也是正態(tài)分布。

用圖來理解這兩條性質(zhì)就是，對于下圖的二維正態(tài)分布 ，單獨的 和 都服從一維正態(tài)分布。

條件概率 的PDF 對應(yīng)圖中的紅線，顯然也是一維正態(tài)分布。

寫一段簡單的代碼驗證二維正態(tài)分布的單個分量服從正態(tài)分布。

代碼中，我們用np.random.normal生成了 10000 個服從二維正態(tài)分布的 x, y 點，然后我們丟棄 y，只保留 x，并畫出 10000 個 x 的分布。

def plot_normal_1d():
    x, _ = np.random.normal(loc=0, scale=1, size=(2, 10000))
    import seaborn as sns
    sns.distplot(x, hist=True, kde=True, bins=100, color='darkblue',
                 hist_kws={'edgecolor': 'black'},
                 kde_kws={'linewidth': 4})
    plt.title('PDF Normal 1D from 2D')
    plt.show()

Box-Muller 原理

雖然無法直接用逆變換方法生成一維正態(tài)分布，但我們卻能通過先生成二維的正態(tài)分布，利用上面一節(jié)的性質(zhì)，生成一維正態(tài)分布。

而 Box-Muller 就是巧妙生成二維正態(tài)分布樣本點的方法。

首先，我們來看看二維正態(tài)分布可以認為是兩個維度是獨立的，每個維度都是正態(tài)分布。此時，其 PDF 可以寫成兩個一維正態(tài)分布 PDF 的乘積。

這種寫法表明，二維正態(tài)分布僅用一個 r 向量就可以充分表達。注意，r 是向量，不僅有大小還有角度，有兩個分量。這兩個分量本質(zhì)上是獨立的，這就是 Box-Muller 方法的巧妙之處。也就是，Box-Muller 通過角度和半徑大小兩個分量的獨立性分別單獨生成并轉(zhuǎn)換成 (x, y) 對。

角度分量是在 范圍均勻采樣，這一點比較直覺好理解。

再來看看半徑分量 r。我們令

則 s 服從指數(shù)分布 。

不信么？我們不妨來做個模擬實驗，下圖是模擬 10000次二維正態(tài)分布 (x, y) 點后轉(zhuǎn)換成 s 的分布。

模擬和plot代碼如下

def plot_r_squared():
    def gen_normal_samples(n):
        x, y = np.random.normal(loc=0, scale=1, size=(2, n))
        return x, y

    x, y = gen_normal_samples(10000)
    s = (x * x + y * y)/2
    plot_dist_1d(s, title='PDF $s = {{x^2 + y^2}\over{2}} \sim exp(1)$')
    

def plot_dist_1d(X, title='PDF '):
    import seaborn as sns
    plt.rcParams.update({
        "text.usetex": True,
        "font.family": "sans-serif",
        "font.sans-serif": ["Helvetica"]})
    sns.distplot(X, hist=True, kde=True, bins=100, color='darkblue',
                 hist_kws={'edgecolor': 'black'},
                 kde_kws={'linewidth': 4})
    plt.title(title)
    plt.show()    

確信了 s 符合指數(shù)分布，根據(jù)指數(shù)分布的 PDF，可以推出二維正態(tài) PDF中的 也符合指數(shù)分布，即

至此，總結(jié)一下Box-Muller方法。我們視二維正態(tài)分布PDF為獨立兩部分的乘積，第一部分是在 范圍中的均勻分布，代表了二維平面中的角度 ，第二部分為 的指數(shù)分布，代表半徑大小。

Box-Muller 方法通過兩個服從 [0, 1] 均勻分布的樣本 u1和u2，轉(zhuǎn)換成獨立的角度和半徑樣本，具體過程如下

1. 生成 [0, 1] 的均勻分布 u1，利用逆變換采樣方法轉(zhuǎn)換成 exp(1) 樣本，此為二維平面點半徑 r

2. 生成 [0, 1] 的均勻分布 u2，乘以 ，即為樣本點的角度

3. 將 r 和  轉(zhuǎn)換成 x, y 坐標(biāo)下的點。

理解了整個過程的意義，下面的代碼就很直白。

def normal_box_muller():
    import random
    from math import sqrt, log, pi, cos, sin
    u1 = random.random()
    u2 = random.random()
    r = sqrt(-2 * log(u1))
    theta = 2 * pi * u2
    z0 = r * cos(theta)
    z1 = r * sin(theta)
    return z0, z1

接下來，我們來看看 Box-Muller 法生成的二維標(biāo)準正態(tài)分布動畫吧。

拒絕采樣極坐標(biāo)方法

Box-Muller 方法還有一種形式，稱為極坐標(biāo)形式，屬于拒絕采樣方法。

1. 生成獨立的 u, v 和 s

分別生成 [0, 1] 均勻分布 u 和 v。令 。如果 s = 0或  s ≥ 1，則丟棄 u 和 v ，并嘗試另一對 (u , v)。因為 u 和 v 是均勻分布的，并且因為只允許單位圓內(nèi)的點，所以 s 的值也將均勻分布在開區(qū)間 (0, 1) 中。注意，這里的 s 的意義雖然也為半徑，但不同于基本方法中的 s。這里 s 取值范圍為 (0, 1) ，目的是通過 s 生成指數(shù)分布，而基本方法中的 s 取值范圍為 [0, +∞]，表示二維正態(tài)分布 PDF 采樣點的半徑。復(fù)用符號 s 的原因是為了對應(yīng)維基百科中關(guān)于基本方法和極坐標(biāo)方法的數(shù)學(xué)描述。

我們用代碼來驗證 s 服從 (0, 1) 范圍上的均勻分布。

def gen_polar_s():
    import random
    while True:
        u = random.uniform(-1, 1)
        v = random.uniform(-1, 1)
        s = u * u + v * v
        if s >= 1.0 or s == 0.0:
            continue
        return s


def plot_polar_s():
    s = [gen_polar_s() for _ in range(10000) ]
    plot_dist_1d(s, title='PDF Polar $s = u^2 + v^2$')

2. 將 u, v, s 轉(zhuǎn)換成 x, y

若將 看成是基本方法中的 u1，就可以用同樣的方式轉(zhuǎn)換成指數(shù)分布，用以代表二維PDF的半徑。

同時，根據(jù)下圖， 和 可以直接用  u, v, R 表示出來，并不需要通過三角函數(shù)顯示計算出 。有了半徑， 和 ，則可以直接計算出 x, y 坐標(biāo)，（下面用 代替 ）。

同樣，Box-Muller 極坐標(biāo)方法的代碼和公式一致。

def normal_box_muller_polar():
    import random
    from math import sqrt, log
    while True:
        u = random.uniform(-1, 1)
        v = random.uniform(-1, 1)
        s = u * u + v * v
        if s >= 1.0 or s == 0.0:
            continue
        z0 = u * sqrt(-2 * log(s) / s)
        z1 = v * sqrt(-2 * log(s) / s)
        return z0, z1

拒絕采樣的效率

極坐標(biāo)方法與基本方法的不同之處在于它是一種拒絕采樣。因此，它會丟棄一些生成的隨機數(shù)，但可能比基本方法更快，因為它計算更簡單：避免使用昂貴的三角函數(shù)，并且在數(shù)值上更穩(wěn)健。極坐標(biāo)方法丟棄了生成總輸入對的 1 ? π /4 ≈ 21.46%，即需要 4/ π ≈ 1.2732 個輸入隨機數(shù)，輸出一個隨機采樣。