服从正态分布python 服从正态分布N(0,1)

服从正态分布随机数的生成

• 生成单变量正态分布随机数

• Box-Muller 算法
• Accept/Reject 算法（接受/拒绝算法）
• 附录

• Box-Muller 算法的简要证明

• 参考

生成单变量正态分布随机数

Box-Muller 算法

Box-Muller算法是利用两个i.i.d. (independent identical distribution)的$U(0, \, 1)$ （(0, 1) 区间的均匀分布）来生成两个i.i.d. 的标准正态分布（$N(0, \, 1)$）的算法。 其具体步骤如下。

1. 首先我们生成两个服从$U(0, \, 1)$的i.i.d.变量，假设为$U_1, \, U_2$。
2. 计算$X_1 = \cos(2\pi U_1) \sqrt{-2 \log(U_2)}$,
$X_2 = \sin(2\pi U_1) \sqrt{-2 \log(U_2)}$
3. 于是$X_1, \, X_2$即为服从$N(0, \, 1)$的两个独立的变量。

我们看到Box-Muller算法实际上是生成了两个服从正态分布的独立变量。当然如果我们只需要生成一个变量，我们取$X_1$即可。利用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 = 10^5$个服从标准正态分布的随机数，其统计分布直方图与理论的pdf (probability density function，即$\displaystyle f(x) = \frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}, \, -\infty < x < \infty$) 比较如下：

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)$非常接近。关于Box-Muller算法的简要证明可见附录。

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

下面我们来看如何用 Accept/Reject 算法来产生服从正态分布的随机数。Accept/Reject算法是一个产生随机数的通用的算法。我们不仅可以用Accept/Reject算法产生正态分布的随机数，也可以产生服从其他分布的随机数。Accept/Reject算法的具体步骤（理论依据）如下 [1]：
假设随机变量 $Y \sim f_Y(y)$, $V \sim f_V(V)$, 即 $Y$的 pdf 为$f_Y(y)$，$V$的 pdf 为$f_V(V)$。$f_Y(y)$ 与 $f_V(V)$ 的定义域相同。并且对任意 $y$，存在 $M$，使得 $\displaystyle \frac{f_Y(y)}{f_V(y)} \leq M < \infty$。于是，为了生成服从 $Y \sim f_Y$的随机变量 $Y$，我们可以采用如下步骤：
a. 生成两个独立的随机变量 $U$ 和 $V$，$U \sim U(0, \, 1)$，即$uniform(0, \, 1)$；$V \sim f_V$。
b. 如果 $U < \frac{1}{M} f_Y(V) / f_V (V)$，我们取 $Y = V$；否则，返回步骤 a。

下面我们证明通过步骤 a 和 b 生成的随机变量 $Y$ 的 pdf 是$f_Y$。

证明的思路是计算概率$P(Y \leq y)$。如果我们能得到$\displaystyle P(Y \leq y) = \int_{-\infty}^y f_Y(y) dy$，那我们就证明了随机变量$Y$的pdf是$f_Y$。
$\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}$ $\square$

于是我们就证明了Accept/Reject 算法。

有了Accept/Reject 算法，我们来看怎么生成标准态分布$N(0, \, 1)$的随机变量。

我们选取 $U \sim U(0, \, 1)$, $V \sim \text{exponential}(\lambda)$。对于指数分布，我们有$f_V(v) = \lambda e^{-\lambda x}, \, x \geq 0$。首先我们须要验证$\displaystyle \frac{f_Y(y)}{f_V(y)} \leq M < \infty, \, \forall y \in \mathbb{R}$ 是否成立。 我们有$\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}}$。于是我们可以取$M$的值为$\displaystyle \frac{1}{\sqrt{2 \pi} \lambda} e^{\frac{\lambda^2}{2}},(\lambda > 0)$。

从以上分析看到对于服从指数分布的随机变量$V \sim \text{exponential}(\lambda)$，我们都可以用$V$与均匀分布$U$相结合来生成服从标准正态分布的随机数。那么，具体取$\lambda$为什么值才能使我们的Accept/Reject算法最高效呢？这里我们取使得$\displaystyle M(\lambda) = \frac{1}{\sqrt{2 \pi} \lambda} e^{\frac{\lambda^2}{2}}$最小的$\lambda$。这是为了使我们的Accept/Reject算法效率最高。即我们能够accept的次数最多 (Reject 的次数最少)。通过对$\displaystyle M(\lambda)$求导我们可知使得$\displaystyle M(\lambda)$最小的$\lambda$为1。从而我们选取$V \sim \text{exponential}(1)$。

有了 $U$ 与 $V$，我们可以根据上述 Accept/Reject 算法的步骤，来生成服从$N(0, \, 1)$的随机变量$Y$，具体步骤如下：
a. 生成独立的两个随机变量$U$和$V$，$U \sim U(0, \, 1)$；$V \sim \text{exponential}(1)$；
b. 如果$U < e^{-\frac{1}{2}(V - 1)^2}$，则取$Y = V$；反之则返回步骤a。注意到指数分布生成的随机变量只能为正，但是我们的正态分布的随机数可正可负。所以我们用一个$(0, \, 1)$均匀分布随机数来决定$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 = 10^5$个服从标准正态分布的随机数，其统计分布直方图与理论的 pdf 比较如下：

附录

Box-Muller 算法的简要证明

对于Box-Muller算法的证明，我们利用两个变量的概率密度函数变换 [1]。我们已知$U_1 \sim U(0, \, 1), \, U_2 \sim U(0, \, 1)$，所以$\displaystyle f_{U_1, \, U_2} (u_1, \, u_2) = 1, \, 0 \leq u_1 \leq 1, \, 0 \leq u_2 \, \leq 1$。我们要求出经过变换$X_1 = \cos(2\pi U_1) \sqrt{-2 \log(U_2)}$,
$X_2 = \sin(2\pi U_1) \sqrt{-2 \log(U_2)}$之后， 新的变量$X_1, \, X_2$的joint distribution。根据两个变量的概率密度函数变换公式，我们有$f_{X_1, \, X_2} (x_1, \, x_2) = 1 \times \vert J \vert$。$J$ 是这个变换的Jacobian，即$\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}$。

为了方便计算，我们用$\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}$。

经过计算，我们有$\displaystyle J = \frac{u_2}{2\pi}$。而把$X_1^2$与$X_2^2$相加，我们有$\displaystyle u_2 =e^{-\frac{x_1^2 + x_2^2}{2}}$ （这里可以看出这个变换是一一对应的）。代入$f_{X_1, \, X_2} (x_1, \, x_2) = 1 \times \vert J \vert$，我们有$\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}}$。从而我们就证明了$X_1, \, X_2$是独立的且服从标准正态分布。

参考

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