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

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

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，即$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$，使得 $\frac{f_Y(y)}{f_V(y)}\le 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\le y)$。如果我们能得到$P(Y\le y)={\int }_{-\infty }^y{f}_Y(y)dy$，那我们就证明了随机变量$Y$的pdf是$f_Y$。
$$\begin{array}{rl}P(Y\le y)& =P\left(V\le y|U<\frac{1}{M}{f}_Y(V)/f_V(V)\right)\\ & =\frac{P\left(V\le 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)}1du}{{\int }_{-\infty }^{\infty }{f}_V(v)dv{\int }_0^{\frac{1}{M}{f}_Y(v)/f_V(v)}1du}\\ & =\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{array}$$ $\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\ge 0$。首先我们须要验证$\frac{f_Y(y)}{f_V(y)}\le M<\infty ,\forall y\in \mathbb{R}$ 是否成立。 我们有$\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}}\le \frac{1}{\sqrt{2\pi }\lambda }{e}^{\frac{{\lambda }^2}{2}}$。于是我们可以取$M$的值为$\frac{1}{\sqrt{2\pi }\lambda }{e}^{\frac{{\lambda }^2}{2}},(\lambda >0)$。

从以上分析看到对于服从指数分布的随机变量$V\sim \text{exponential}(\lambda )$，我们都可以用$V$与均匀分布$U$相结合来生成服从标准正态分布的随机数。那么，具体取$\lambda$为什么值才能使我们的Accept/Reject算法最高效呢？这里我们取使得$M(\lambda )=\frac{1}{\sqrt{2\pi }\lambda }{e}^{\frac{{\lambda }^2}{2}}$最小的$\lambda$。这是为了使我们的Accept/Reject算法效率最高。即我们能够accept的次数最多 (Reject 的次数最少)。通过对$M(\lambda )$求导我们可知使得$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)$，所以$f_{U_1,U_2}(u_1,u_2)=1,0\le u_1\le 1,0\le u_2\le 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 |J|$。$J$ 是这个变换的Jacobian，即$J=\left|\begin{array}{cc}\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{array}\right|$。

为了方便计算，我们用$J=1/\left|\begin{array}{cc}\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{array}\right|$。

经过计算，我们有$J=\frac{u_2}{2\pi }$。而把$X_1^2$与$X_2^2$相加，我们有$u_2=e^{-\frac{x_1^2+x_2^2}{2}}$ （这里可以看出这个变换是一一对应的）。代入$f_{X_1,X_2}(x_1,x_2)=1\times |J|$，我们有$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