服從正态分布随機數的生成

• 生成單變量正态分布随機數
• • Box-Muller 算法
  • Accept/Reject 算法（接受/拒絕算法）
  • 附錄
  • • Box-Muller 算法的簡要證明
  • 參考

生成單變量正态分布随機數

Box-Muller 算法

Box-Muller算法是利用兩個i.i.d. (independent identical distribution)的 U ( 0 ,   1 ) U(0, \, 1) U(0,1) （(0, 1) 區間的均勻分布）來生成兩個i.i.d. 的标準正态分布（ N ( 0 ,   1 ) N(0, \, 1) N(0,1)）的算法。 其具體步驟如下。

1. 首先我們生成兩個服從 U ( 0 ,   1 ) U(0, \, 1) U(0,1)的i.i.d.變量，假設為 U 1 ,   U 2 U_1, \, U_2 U1​,U2​。

2. 計算 X 1 = cos ⁡ ( 2 π U 1 ) − 2 log ⁡ ( U 2 ) X_1 = \cos(2\pi U_1) \sqrt{-2 \log(U_2)} X1​=cos(2πU1​)−2log(U2​)

​,

X 2 = sin ⁡ ( 2 π U 1 ) − 2 log ⁡ ( U 2 ) X_2 = \sin(2\pi U_1) \sqrt{-2 \log(U_2)} X2​=sin(2πU1​)−2log(U2​)

​

3. 于是 X 1 ,   X 2 X_1, \, X_2 X1​,X2​即為服從 N ( 0 ,   1 ) N(0, \, 1) N(0,1)的兩個獨立的變量。

我們看到Box-Muller算法實際上是生成了兩個服從正态分布的獨立變量。當然如果我們隻需要生成一個變量，我們取 X 1 X_1 X1​即可。利用Python，Box-Muller算法的代碼如下：

import numpy as np
from scipy.stats import norm
from scipy.stats import uniform
import matplotlib.pyplot as plt
import pandas as pd

class Box_Muller:
    
    def __init__(self, n: int):
        """
        Suppose we want to generate n i.i.d. N(0, 1) random variables
        """
        self.N = n
    
    def generate_Box_Muller(self) -> 'list(float)':
        """
        Use Box-Muller algorithm to generate self.N number of i.i.d. N(0, 1) random variables.
        """
        res = []
        for i in range(n):
            U1, U2 = np.random.uniform(), np.random.uniform()
            X1, X2 = np.cos(2 * np.pi * U1) * np.sqrt(-2 * np.log(U2)), \
                    np.sin(2 * np.pi * U1) * np.sqrt(-2 * np.log(U2))
            res.append(X1)
        return res
           

如果我們生成 n = 1 0 5 n = 10^5 n=105個服從标準正态分布的随機數，其統計分布直方圖與理論的pdf (probability density function，即 f ( x ) = 1 2 π e − x 2 2 ,   − ∞ < x < ∞ \displaystyle f(x) = \frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}, \, -\infty < x < \infty f(x)=2π

​1​e2−x2​,−∞<x<∞) 比較如下：

n = 10 ** 5
a = Box_Muller(n)
x = a.generate_Box_Muller()

plt.figure(figsize=(8, 6), dpi=100)
plt.hist(x, bins = 100, density=True)
x_forplot = np.linspace(-5, 5, 1000)
plt.plot(x_forplot, norm.pdf(x_forplot), linewidth = 3)
plt.xlabel("random variable values", fontsize=20)
plt.ylabel("histogram frequency", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.legend(['theoretical pdf', 'Bos-Muller sample frequencies'], fontsize=10)
           

我們可以看到統計分布直方圖與理論 f ( x ) f(x) f(x)非常接近。關于Box-Muller算法的簡要證明可見附錄。

Accept/Reject 算法（接受/拒絕算法）

下面我們來看如何用 Accept/Reject 算法來産生服從正态分布的随機數。Accept/Reject算法是一個産生随機數的通用的算法。我們不僅可以用Accept/Reject算法産生正态分布的随機數，也可以産生服從其他分布的随機數。Accept/Reject算法的具體步驟（理論依據）如下 [1]：

假設随機變量 Y ∼ f Y ( y ) Y \sim f_Y(y) Y∼fY​(y), V ∼ f V ( V ) V \sim f_V(V) V∼fV​(V), 即 Y Y Y的 pdf 為 f Y ( y ) f_Y(y) fY​(y)， V V V的 pdf 為 f V ( V ) f_V(V) fV​(V)。 f Y ( y ) f_Y(y) fY​(y) 與 f V ( V ) f_V(V) fV​(V) 的定義域相同。并且對任意 y y y，存在 M M M，使得 f Y ( y ) f V ( y ) ≤ M < ∞ \displaystyle \frac{f_Y(y)}{f_V(y)} \leq M < \infty fV​(y)fY​(y)​≤M<∞。于是，為了生成服從 Y ∼ f Y Y \sim f_Y Y∼fY​的随機變量 Y Y Y，我們可以采用如下步驟：

a. 生成兩個獨立的随機變量 U U U 和 V V V， U ∼ U ( 0 ,   1 ) U \sim U(0, \, 1) U∼U(0,1)，即 u n i f o r m ( 0 ,   1 ) uniform(0, \, 1) uniform(0,1)； V ∼ f V V \sim f_V V∼fV​。

b. 如果 U < 1 M f Y ( V ) / f V ( V ) U < \frac{1}{M} f_Y(V) / f_V (V) U<M1​fY​(V)/fV​(V)，我們取 Y = V Y = V Y=V；否則，傳回步驟 a。

下面我們證明通過步驟 a 和 b 生成的随機變量 Y Y Y 的 pdf 是 f Y f_Y fY​。

證明的思路是計算機率 P ( Y ≤ y ) P(Y \leq y) P(Y≤y)。如果我們能得到 P ( Y ≤ y ) = ∫ − ∞ y f Y ( y ) d y \displaystyle P(Y \leq y) = \int_{-\infty}^y f_Y(y) dy P(Y≤y)=∫−∞y​fY​(y)dy，那我們就證明了随機變量 Y Y Y的pdf是 f Y f_Y fY​。

P ( Y ≤ y ) = P ( V ≤ y   ∣   U < 1 M f Y ( V ) / f V ( V ) ) = P ( V ≤ y ,   U < 1 M f Y ( V ) / f V ( V ) ) P ( U < 1 M f Y ( V ) / f V ( V ) ) = ∫ − ∞ y f V ( v ) d v ∫ 0 1 M f Y ( v ) / f V ( v ) 1   d u ∫ − ∞ ∞ f V ( v ) d v ∫ 0 1 M f Y ( v ) / f V ( v ) 1   d u = ∫ − ∞ y f V ( v ) d v 1 M f Y ( v ) / f V ( v ) ∫ − ∞ ∞ f V ( v ) d v 1 M f Y ( v ) / f V ( v ) = 1 M ∫ − ∞ y f Y ( v ) d v 1 M ∫ − ∞ ∞ f Y ( v ) d v = ∫ − ∞ y f Y ( y ) d y . ( note that ∫ − ∞ ∞ f Y ( v ) d v = 1 ) \begin{aligned} P(Y \leq y) &= P \left(V \leq y \, \vert \, U < \frac{1}{M} f_Y(V) / f_V (V) \right) \\ &= \frac{P \left(V \leq y, \, U < \frac{1}{M} f_Y(V) / f_V (V) \right) }{P \left( U < \frac{1}{M} f_Y(V) / f_V (V) \right)} \\ &= \frac{\int_{-\infty}^{y} f_V(v) dv \int_{0}^{\frac{1}{M} f_Y(v) / f_V (v)} 1 \,du}{\int_{-\infty}^{\infty} f_V(v) dv \int_{0}^{\frac{1}{M} f_Y(v) / f_V (v)} 1 \,du} \\ &= \frac{\int_{-\infty}^{y} f_V(v) dv \frac{1}{M} f_Y(v) / f_V (v)}{\int_{-\infty}^{\infty} f_V(v) dv \frac{1}{M} f_Y(v) / f_V (v)} \\ &= \frac{\frac{1}{M} \int_{-\infty}^{y} f_Y(v) dv}{\frac{1}{M} \int_{-\infty}^{\infty} f_Y(v) dv} \\ &= \int_{-\infty}^y f_Y(y) dy. \\ & (\text{note that} \int_{-\infty}^{\infty} f_Y(v) dv = 1) \end{aligned} P(Y≤y)​=P(V≤y∣U<M1​fY​(V)/fV​(V))=P(U<M1​fY​(V)/fV​(V))P(V≤y,U<M1​fY​(V)/fV​(V))​=∫−∞∞​fV​(v)dv∫0M1​fY​(v)/fV​(v)​1du∫−∞y​fV​(v)dv∫0M1​fY​(v)/fV​(v)​1du​=∫−∞∞​fV​(v)dvM1​fY​(v)/fV​(v)∫−∞y​fV​(v)dvM1​fY​(v)/fV​(v)​=M1​∫−∞∞​fY​(v)dvM1​∫−∞y​fY​(v)dv​=∫−∞y​fY​(y)dy.(note that∫−∞∞​fY​(v)dv=1)​ □ \square □

于是我們就證明了Accept/Reject 算法。

有了Accept/Reject 算法，我們來看怎麼生成标準态分布 N ( 0 ,   1 ) N(0, \, 1) N(0,1)的随機變量。

我們選取 U ∼ U ( 0 ,   1 ) U \sim U(0, \, 1) U∼U(0,1), V ∼ exponential ( λ ) V \sim \text{exponential}(\lambda) V∼exponential(λ)。對于指數分布，我們有 f V ( v ) = λ e − λ x ,   x ≥ 0 f_V(v) = \lambda e^{-\lambda x}, \, x \geq 0 fV​(v)=λe−λx,x≥0。首先我們須要驗證 f Y ( y ) f V ( y ) ≤ M < ∞ ,   ∀ y ∈ R \displaystyle \frac{f_Y(y)}{f_V(y)} \leq M < \infty, \, \forall y \in \mathbb{R} fV​(y)fY​(y)​≤M<∞,∀y∈R 是否成立。 我們有 f Y ( y ) f V ( y ) = 1 2 π e − x 2 2 λ e − λ x = 1 2 π λ e − 1 2 ( y − λ ) 2 + λ 2 2 ≤ 1 2 π λ e λ 2 2 \displaystyle \frac{f_Y(y)}{f_V(y)} = \frac{\frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}}{ \lambda e^{-\lambda x}} = \frac{1}{\sqrt{2 \pi} \lambda} e^{-\frac{1}{2}(y - \lambda)^2 + \frac{\lambda^2}{2}} \leq \frac{1}{\sqrt{2 \pi} \lambda} e^{\frac{\lambda^2}{2}} fV​(y)fY​(y)​=λe−λx2π

​1​e2−x2​​=2π

​λ1​e−21​(y−λ)2+2λ2​≤2π

​λ1​e2λ2​。于是我們可以取 M M M的值為 1 2 π λ e λ 2 2 , ( λ > 0 ) \displaystyle \frac{1}{\sqrt{2 \pi} \lambda} e^{\frac{\lambda^2}{2}},(\lambda > 0) 2π

​λ1​e2λ2​,(λ>0)。

從以上分析看到對于服從指數分布的随機變量 V ∼ exponential ( λ ) V \sim \text{exponential}(\lambda) V∼exponential(λ)，我們都可以用 V V V與均勻分布 U U U相結合來生成服從标準正态分布的随機數。那麼，具體取 λ \lambda λ為什麼值才能使我們的Accept/Reject算法最高效呢？這裡我們取使得 M ( λ ) = 1 2 π λ e λ 2 2 \displaystyle M(\lambda) = \frac{1}{\sqrt{2 \pi} \lambda} e^{\frac{\lambda^2}{2}} M(λ)=2π

​λ1​e2λ2​最小的 λ \lambda λ。這是為了使我們的Accept/Reject算法效率最高。即我們能夠accept的次數最多 (Reject 的次數最少)。通過對 M ( λ ) \displaystyle M(\lambda) M(λ)求導我們可知使得 M ( λ ) \displaystyle M(\lambda) M(λ)最小的 λ \lambda λ為1。進而我們選取 V ∼ exponential ( 1 ) V \sim \text{exponential}(1) V∼exponential(1)。

有了 U U U 與 V V V，我們可以根據上述 Accept/Reject 算法的步驟，來生成服從 N ( 0 ,   1 ) N(0, \, 1) N(0,1)的随機變量 Y Y Y，具體步驟如下：

a. 生成獨立的兩個随機變量 U U U和 V V V， U ∼ U ( 0 ,   1 ) U \sim U(0, \, 1) U∼U(0,1)； V ∼ exponential ( 1 ) V \sim \text{exponential}(1) V∼exponential(1)；

b. 如果 U < e − 1 2 ( V − 1 ) 2 U < e^{-\frac{1}{2}(V - 1)^2} U<e−21​(V−1)2，則取 Y = V Y = V Y=V；反之則傳回步驟a。注意到指數分布生成的随機變量隻能為正，但是我們的正态分布的随機數可正可負。是以我們用一個 ( 0 ,   1 ) (0, \, 1) (0,1)均勻分布随機數來決定 V V V的符号。即如果随機數小于0.5，則取為正；反之則取為負。

根據以上Accept/Reject算法的Python代碼如下：

class Accept_Reject:
    
    def __init__(self, n: int):
        """
        Suppose we want to generate n i.i.d. N(0, 1) random variables
        """
        self.N = n
    
    def generate_accept_reject(self) -> 'list(float)':
        """
        Use accept/reject algorithm to generate self.N random numbers that follows 
        N(0, 1) distribution.
        """
        standNorm = []
        M = 1 / np.sqrt(np.pi * 2) * np.e ** (0.5)
        while len(standNorm) < self.N:
            U = np.random.uniform()
            W = np.random.uniform()
            V = -np.log(W)
            if U < np.e ** (-0.5 * (V - 1) ** 2):
                if np.random.uniform() <= 0.5:
                    standNorm.append(V)
                else:
                    standNorm.append(-V)
        return standNorm
           

如果我們用Accept/Reject生成 n = 1 0 5 n = 10^5 n=105個服從标準正态分布的随機數，其統計分布直方圖與理論的 pdf 比較如下：

附錄

Box-Muller 算法的簡要證明

對于Box-Muller算法的證明，我們利用兩個變量的機率密度函數變換 [1]。我們已知 U 1 ∼ U ( 0 ,   1 ) ,   U 2 ∼ U ( 0 ,   1 ) U_1 \sim U(0, \, 1), \, U_2 \sim U(0, \, 1) U1​∼U(0,1),U2​∼U(0,1)，是以 f U 1 ,   U 2 ( u 1 ,   u 2 ) = 1 ,   0 ≤ u 1 ≤ 1 ,   0 ≤ u 2   ≤ 1 \displaystyle f_{U_1, \, U_2} (u_1, \, u_2) = 1, \, 0 \leq u_1 \leq 1, \, 0 \leq u_2 \, \leq 1 fU1​,U2​​(u1​,u2​)=1,0≤u1​≤1,0≤u2​≤1。我們要求出經過變換 X 1 = cos ⁡ ( 2 π U 1 ) − 2 log ⁡ ( U 2 ) X_1 = \cos(2\pi U_1) \sqrt{-2 \log(U_2)} X1​=cos(2πU1​)−2log(U2​)

​,

X 2 = sin ⁡ ( 2 π U 1 ) − 2 log ⁡ ( U 2 ) X_2 = \sin(2\pi U_1) \sqrt{-2 \log(U_2)} X2​=sin(2πU1​)−2log(U2​)

​之後， 新的變量 X 1 ,   X 2 X_1, \, X_2 X1​,X2​的joint distribution。根據兩個變量的機率密度函數變換公式，我們有 f X 1 ,   X 2 ( x 1 ,   x 2 ) = 1 × ∣ J ∣ f_{X_1, \, X_2} (x_1, \, x_2) = 1 \times \vert J \vert fX1​,X2​​(x1​,x2​)=1×∣J∣。 J J J 是這個變換的Jacobian，即 J = ∣ ∂ u 1 ∂ x 1 ∂ u 1 ∂ x 2 ∂ u 2 ∂ x 1 ∂ u 2 ∂ x 2 ∣ \displaystyle J = \begin{vmatrix} \frac{\partial u_1}{\partial x_1} & \frac{\partial u_1}{\partial x_2} \\ \\ \frac{\partial u_2}{\partial x_1} & \frac{\partial u_2}{\partial x_2} \end{vmatrix} J=∣∣∣∣∣∣​∂x1​∂u1​​∂x1​∂u2​​​∂x2​∂u1​​∂x2​∂u2​​​∣∣∣∣∣∣​。

為了友善計算，我們用 J = 1 / ∣ ∂ x 1 ∂ u 1 ∂ x 1 ∂ u 2 ∂ x 2 ∂ u 1 ∂ x 2 ∂ u 2 ∣ \displaystyle J = 1 / \begin{vmatrix} \frac{\partial x_1}{\partial u_1} & \frac{\partial x_1}{\partial u_2} \\ \frac{\partial x_2}{\partial u_1} & \frac{\partial x_2}{\partial u_2} \end{vmatrix} J=1/∣∣∣∣∣​∂u1​∂x1​​∂u1​∂x2​​​∂u2​∂x1​​∂u2​∂x2​​​∣∣∣∣∣​。

經過計算，我們有 J = u 2 2 π \displaystyle J = \frac{u_2}{2\pi} J=2πu2​​。而把 X 1 2 X_1^2 X12​與 X 2 2 X_2^2 X22​相加，我們有 u 2 = e − x 1 2 + x 2 2 2 \displaystyle u_2 =e^{-\frac{x_1^2 + x_2^2}{2}} u2​=e−2x12​+x22​​ （這裡可以看出這個變換是一一對應的）。代入 f X 1 ,   X 2 ( x 1 ,   x 2 ) = 1 × ∣ J ∣ f_{X_1, \, X_2} (x_1, \, x_2) = 1 \times \vert J \vert fX1​,X2​​(x1​,x2​)=1×∣J∣，我們有 f X 1 ,   X 2 ( x 1 ,   x 2 ) = 1 2 π e − x 1 2 + x 2 2 2 = 1 2 π e − x 1 2 2 × 1 2 π e − x 2 2 2 \displaystyle f_{X_1, \, X_2} (x_1, \, x_2) = \frac{1}{2\pi} e^{-\frac{x_1^2 + x_2^2}{2}} = \frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}} \times \frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}} fX1​,X2​​(x1​,x2​)=2π1​e−2x12​+x22​​=2π

​1​e−2x12​​×2π

​1​e−2x22​​。進而我們就證明了 X 1 ,   X 2 X_1, \, X_2 X1​,X2​是獨立的且服從标準正态分布。

參考

[1] George Casella, Roger L. Berger, Statistical inference, Chapter 4.3