Box-Muller 算法可以把独立的均匀分布随机数变换为独立的正态分布随机数

设 U1,U2∼i.i.d.U(0,1)U_1, U_2 \stackrel{\text{i.i.d.}}{\sim} U(0,1)，则 −2ln⁡U1cos⁡(2πU2)\sqrt{-2\ln U_1} \cos(2\pi U_2) 和 −2ln⁡U1sin⁡(2πU2)\sqrt{-2\ln U_1} \sin(2\pi U_2) 独立且都服从标准正态分布 N(0,1)\mathcal{N}(0,1)。

这种算法简单得令人难以置信，不过证明它的正确性其实并不困难：核心思路是对二维正态分布进行极坐标变换。只要注意到二维正态分布独特的对称性质，那么接下来的一切推导就都是普通和理所当然的。

数学证明

二维正态分布的极坐标表示

众所周知，一维正态分布的概率密度函数为

fX(x)=12πσe−(x−μ)22σ2x∈R f_X(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{\frac{-(x-\mu)^2}{2\sigma^2}} \quad x\in\mathbb{R}

令 μ=0,σ=1\mu=0, \sigma=1 就可以得到标准正态分布

fX(x)=12πe−x22 f_X(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}

二维正态分布与一维的类似。设 X,Y∼i.i.d.N(0,1)X, Y \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0,1)，则其联合概率密度函数为

fX,Y(x,y)=12πe−x2+y22(x,y)∈R2 f_{X,Y}(x,y) = \frac{1}{2\pi} e^{-\frac{x^2+y^2}{2}} \quad (x,y)\in\mathbb{R}^2

对随机变量作极坐标变换

{X=Rcos⁡ΘY=Rsin⁡Θ \begin{cases} X = R\cos\Theta \\ Y = R\sin\Theta \end{cases}

得到

{x(r,θ)=rcos⁡θy(r,θ)=rsin⁡θ \begin{cases} x(r, \theta) = r\cos\theta \\ y(r, \theta) = r\sin\theta \end{cases}

其中

{r≥0θ∈[0,2π) \begin{cases} r\ge 0 \\ \theta\in[0,2\pi) \end{cases}

变换的雅可比矩阵为

J=[∂x∂r∂x∂θ∂y∂r∂y∂θ]=[cos⁡θ−rsin⁡θsin⁡θrcos⁡θ] J = \begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \end{bmatrix} = \begin{bmatrix} \cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta \end{bmatrix}

雅可比行列式为

∣J∣=∣cos⁡θ−rsin⁡θsin⁡θrcos⁡θ∣=rcos⁡2θ+rsin⁡2θ=r |J| = \begin{vmatrix} \cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta \end{vmatrix} = r\cos^2\theta + r\sin^2\theta = r

得到 R,ΘR,\Theta 的联合概率密度

fR,Θ(r,θ)=fX,Y(x(r,θ),y(r,θ))∣J∣=12πe−(rcos⁡θ)2+(rsin⁡θ)22r=12πre−r22 \begin{align*} f_{R,\Theta}(r,\theta) &= f_{X,Y}(x(r,\theta),y(r, \theta)) |J| \\ &= \frac{1}{2\pi} e^{-\frac{(r\cos\theta)^2+(r\sin\theta)^2}{2}}r \\ &= \frac{1}{2\pi} r e^{-\frac{r^2}{2}} \end{align*}

即

fR,Θ(r,θ)=12πre−r22r≥0,θ∈[0,2π) f_{R,\Theta}(r,\theta) = \frac{1}{2\pi} r e^{-\frac{r^2}{2}} \quad r\ge 0, \theta\in[0,2\pi)

计算 Θ\Theta 的边缘密度

fΘ(θ)=∫0∞fR,Θ(r,θ)dr=∫0∞12πre−r22dr f_\Theta(\theta) = \int_0^{\infty} f_{R,\Theta}(r,\theta) \mathrm{d}r = \int_0^{\infty} \frac{1}{2\pi} r e^{-\frac{r^2}{2}} \mathrm{d}r

不过这里并不需要真的进行积分——显然积分的结果是个和 θ\theta 无关的常数。于是 fΘ(θ)f_\Theta(\theta) 在取值范围内为常数。又已知 θ∈[0,2π)\theta\in[0,2\pi)，所以

Θ∼U(0,2π) \Theta \sim U(0,2\pi)

计算 RR 的边缘密度

fR(r)=∫02πfR,Θ(r,θ)dθ=∫02π12πre−r22dθ=re−r22 f_R(r) = \int_0^{2\pi} f_{R,\Theta}(r,\theta) \mathrm{d}\theta = \int_0^{2\pi} \frac{1}{2\pi} r e^{-\frac{r^2}{2}} \mathrm{d}\theta = r e^{-\frac{r^2}{2}}

如果缺乏注意力的话，这里不太能够直接看出 RR 的分布。但这个概率密度函数和指数分布的 fX(x)=λe−λxf_X(x)=\lambda e^{-\lambda x} 有点相似，因此接下来再对 RR 进行换元。

令 R=SR = \sqrt{S}，则 r(s)=sr(s) = \sqrt{s}，从而有

fS(s)=fR(r(s))∣drds∣=fR(s)∣dsds∣=se−s212s=12e−s2 \begin{align*} f_S(s) &= f_R(r(s)) \left| \frac{\mathrm{d}r}{\mathrm{d}s} \right| \\ &= f_R(\sqrt{s}) \left| \frac{\mathrm{d}\sqrt{s}}{\mathrm{d}s} \right| \\ &= \sqrt{s} e^{-\frac{s}{2}} \frac{1}{2\sqrt{s}} \\ &= \frac{1}{2} e^{-\frac{s}{2}} \end{align*}

现在和指数分布的概率密度完全一样了。因此 SS 服从 λ=1/2\lambda=1/2 的指数分布，即

S∼Exp⁡(1/2) S \sim \operatorname{Exp}(1/2)

下一步就是把均匀分布的随机变量变为指数分布的随机变量。

随机变量逆变换

将随机变量 XX 变为 YY 的通用公式为

Y=FY−1(FX(X)) Y=F_Y^{-1}(F_X(X))

该公式分为两个部分。首先对源分布 XX 进行概率积分变换 FX(X)F_X(X) 得到均匀分布 UU

P(FX(X)≤u)=P(X≤FX−1(u))=FX(FX−1(u))=u=P(U≤u) P(\textcolor{red}{F_X(X)} \le u) = P(X \le F_X^{-1}(u)) = F_X(F_X^{-1}(u)) = u = P(\textcolor{red}{U} \le u)

再对均匀分布 UU 进行逆变换采样 FY−1(U)F_Y^{-1}(U) 得到目标分布 YY

P(FY−1(U)≤y)=P(U≤FY(y))=FY(y)=P(Y≤y) P(\textcolor{red}{F_Y^{-1}(U)} \le y) = P(U \le F_Y(y)) = F_Y(y) = P(\textcolor{red}{Y} \le y)

由于源分布已经是均匀分布 U(0,1)U(0,1)，因此接下来只需计算指数分布的累积分布函数 FSF_S ，然后再求其反函数 FS−1F_S^{-1} 即可。

对于 Exp⁡(1/2)\operatorname{Exp}(1/2)，其累积分布函数为

FS(s)=∫0sfS(s)ds=∫0s12e−s2ds=[−e−t2]0s=1−e−s2 \begin{align*} F_S(s) &= \int_0^{s} f_S(s) \mathrm{d}s \\ &= \int_0^{s} \frac{1}{2} e^{-\frac{s}{2}} \mathrm{d}s \\ &= [-e^{-\frac{t}{2}}]_0^s \\ &= 1 - e^{-\frac{s}{2}} \end{align*}

推导其反函数

e−s2=1−FS−s2=ln⁡(1−FS)s=−2ln⁡(1−FS) \begin{align*} e^{-\frac{s}{2}} &= 1 - F_S \\ -\frac{s}{2} &= \ln(1 - F_S) \\ s &= -2\ln(1 - F_S) \end{align*}

即

FS−1(u)=−2ln⁡(1−u)u∈[0,1) F_S^{-1}(u)=-2\ln(1-u) \quad u\in[0,1)

设 U1∼U(0,1)U_1 \sim U(0,1)，于是有

S=FS−1(U1)=−2ln⁡(1−U1) S = F_S^{-1}(U_1)=-2\ln(1-U_1)

由 1−U1∼U(0,1)1-U_1 \sim U(0,1) 进一步得到

S=−2ln⁡U1 S = -2\ln U_1

而之前已经得到 Θ∼U(0,2π)\Theta \sim U(0,2\pi)。设 U2∼U(0,1)U_2 \sim U(0,1)，则

Θ=2πU2 \Theta = 2\pi U_2

于是

{R=S=−2ln⁡U1Θ=2πU2 \begin{cases} R = \sqrt{S} = \sqrt{-2\ln U_1} \\ \Theta = 2\pi U_2 \end{cases}

带入

{X=Rcos⁡ΘY=Rsin⁡Θ \begin{cases} X = R\cos\Theta \\ Y = R\sin\Theta \end{cases}

得到最终结果

X=−2ln⁡U1cos⁡(2πU2)Y=−2ln⁡U1sin⁡(2πU2) \begin{align*} X &= \sqrt{-2\ln U_1}\cos(2\pi U_2) \\ Y &= \sqrt{-2\ln U_1}\sin(2\pi U_2) \end{align*}

显然，若 U1U_1 与 U2U_2 独立，则 XX 与 YY 独立。这样就成功用一对独立的均匀分布随机数生成了独立的高斯分布随机数。

代码实现

由于我就是在 最近做的一个 Web 项目 中用到了 Box-Muller 算法（这也是写本文的原因），因此这里给出该算法的 JavaScript 实现。

其实 JavaScript 的 Math.random() 能够生成 U(0,1)U(0,1) 随机数，但该方法不允许设置种子，导致缺少很多乐趣，所以我这里先实现了一个基于线性同余的均匀随机数生成器。

function createUniformRandomGenerator(seed) {
  let state = seed >>> 0;

  return () => {
    state = (1664525 * state + 1013904223) >>> 0;
    return state / 0xffffffff;
  };
}

function createGaussianRandomGenerator(seed) {
  const getUniformRandom = createUniformRandomGenerator(seed);
  let spare = null;

  return () => {
    if (spare !== null) {
      const value = spare;
      spare = null;
      return value;
    }

    let u_1 = getUniformRandom();
    while (u_1 === 0) {
      u_1 = getUniformRandom();
    }
    const u_2 = getUniformRandom();

    const r = Math.sqrt(-2 * Math.log(u_1));
    const theta = 2 * Math.PI * u_2;

    spare = r * Math.sin(theta);
    return r * Math.cos(theta);
  };
}

实现中有两个需要注意的点

  1. 由于不能对 00 取对数，因此当 u_1 随机到 00 时需要重试
  2. 由于算法一次调用生成两个随机数，但只返回一个，所以用 spare 缓存另一个，这样偶数次调用时就不会重复计算了

调用示例

SEED=114514
COUNT=5

const g = createGaussianRandomGenerator(SEED)
for (let i = 0; i < COUNT; i++) console.log(g())

输出为

0.6756921337282962
0.7153015025880592
-1.090593327207687
-0.43173865912412884
2.1300190649693405

本文选择给出 JavaScript 实现，其实还有一个好处：每个能浏览网页的设备，事实上都有 JavaScript 的运行环境。对于桌面端，只需在浏览器里按下 F12 打开 DevTools，然后找到 控制台（Console）面板，再把上面的两段代码复制进去就能运行了

浏览器控制台截图

移动端的话没法打开 DevTools，可以到能够在线运行 JavaScript 的网站上尝试。