【數(shù)學(xué)基礎(chǔ)】詳解Box-Muller方法生成正態(tài)分布

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

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

嘗試逆變換方法

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

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

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

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

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

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

二維映射到一維

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

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

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

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

代碼中，我們用np.random.normal生成了 10000 個(gè)服從二維正態(tài)分布的 x, y 點(diǎn)，然后我們丟棄 y，只保留 x，并畫出 10000 個(gè) 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 原理

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

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

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

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

角度分量是在 范圍均勻采樣，這一點(diǎn)比較直覺(jué)好理解。

再來(lái)看看半徑分量 r。我們令

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

不信么？我們不妨來(lái)做個(gè)模擬實(shí)驗(yàn)，下圖是模擬 10000次二維正態(tài)分布 (x, y) 點(diǎn)后轉(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為獨(dú)立兩部分的乘積，第一部分是在 范圍中的均勻分布，代表了二維平面中的角度 ，第二部分為 的指數(shù)分布，代表半徑大小。

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

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

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

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

理解了整個(gè)過(guò)程的意義，下面的代碼就很直白。

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

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

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

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

1. 生成獨(dú)立的 u, v 和 s

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

我們用代碼來(lái)驗(yàn)證 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的半徑。

同時(shí)，根據(jù)下圖， 和 可以直接用  u, v, R 表示出來(lái)，并不需要通過(guò)三角函數(shù)顯示計(jì)算出 。有了半徑， 和 ，則可以直接計(jì)算出 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)方法與基本方法的不同之處在于它是一種拒絕采樣。因此，它會(huì)丟棄一些生成的隨機(jī)數(shù)，但可能比基本方法更快，因?yàn)樗?jì)算更簡(jiǎn)單：避免使用昂貴的三角函數(shù)，并且在數(shù)值上更穩(wěn)健。極坐標(biāo)方法丟棄了生成總輸入對(duì)的 1 ? π /4 ≈ 21.46%，即需要 4/ π ≈ 1.2732 個(gè)輸入隨機(jī)數(shù)，輸出一個(gè)隨機(jī)采樣。