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

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

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π 1e2x2,<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)

Box-Muller 比较
我们可以看到统计分布直方图与理论 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) YfY(y), V ∼ f V ( V ) V \sim f_V(V) VfV(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 YfY的随机变量 Y Y Y,我们可以采用如下步骤:
a. 生成两个独立的随机变量 U U U V V V U ∼ U ( 0 ,   1 ) U \sim U(0, \, 1) UU(0,1),即 u n i f o r m ( 0 ,   1 ) uniform(0, \, 1) uniform(0,1) V ∼ f V V \sim f_V VfV
b. 如果 U < 1 M f Y ( V ) / f V ( V ) U < \frac{1}{M} f_Y(V) / f_V (V) U<M1fY(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(Yy)。如果我们能得到 P ( Y ≤ y ) = ∫ − ∞ y f Y ( y ) d y \displaystyle P(Y \leq y) = \int_{-\infty}^y f_Y(y) dy P(Yy)=yfY(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(Yy)=P(VyU<M1fY(V)/fV(V))=P(U<M1fY(V)/fV(V))P(Vy,U<M1fY(V)/fV(V))=fV(v)dv0M1fY(v)/fV(v)1duyfV(v)dv0M1fY(v)/fV(v)1du=fV(v)dvM1fY(v)/fV(v)yfV(v)dvM1fY(v)/fV(v)=M1fY(v)dvM1yfY(v)dv=yfY(y)dy.(note thatfY(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) UU(0,1), V ∼ exponential ( λ ) V \sim \text{exponential}(\lambda) Vexponential(λ)。对于指数分布,我们有 f V ( v ) = λ e − λ x ,   x ≥ 0 f_V(v) = \lambda e^{-\lambda x}, \, x \geq 0 fV(v)=λeλx,x0。首先我们须要验证 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<,yR 是否成立。 我们有 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π 1e2x2=2π λ1e21(yλ)2+2λ22π λ1e2λ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π λ1e2λ2,(λ>0)

从以上分析看到对于服从指数分布的随机变量 V ∼ exponential ( λ ) V \sim \text{exponential}(\lambda) Vexponential(λ),我们都可以用 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π λ1e2λ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) Vexponential(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) UU(0,1) V ∼ exponential ( 1 ) V \sim \text{exponential}(1) Vexponential(1)
b. 如果 U < e − 1 2 ( V − 1 ) 2 U < e^{-\frac{1}{2}(V - 1)^2} U<e21(V1)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 比较如下:

Accept_Reject

附录

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) U1U(0,1),U2U(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,0u11,0u21。我们要求出经过变换 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=x1u1x1u2x2u1x2u2

为了方便计算,我们用 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/u1x1u1x2u2x1u2x2

经过计算,我们有 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=e2x12+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π1e2x12+x22=2π 1e2x12×2π 1e2x22。从而我们就证明了 X 1 ,   X 2 X_1, \, X_2 X1,X2是独立的且服从标准正态分布。

参考

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

评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值