Box-Muller 变换：从均匀分布生成正态分布的推导过程

准备：你手头有均匀分布随机数生成器，能产生 U(0,1) 的独立样本。目标是用这些均匀样本生成服从标准正态分布 N(0,1) 的独立随机数。

1. 理解逆变换法的局限

回忆：若已知累积分布函数 F(x) 及其反函数 F^{-1}，则令 U ~ U(0,1)，X = F^{-1}(U) 服从 F 对应的分布。但正态分布的 F^{-1} 没有封闭形式（需用 erf 函数），直接计算昂贵。

引入：Box-Muller 变换绕过逆变换，通过两个均匀数生成两个独立的正态数，速度更快。

2. 定义两个独立的均匀随机数

生成：两个独立同分布的均匀随机变量 U_1 和 U_2，均服从 U(0,1)。

构造：利用极坐标思想，定义

• 半径平方：R = sqrt(-2 ln U_1) （注意保证 R ≥ 0）
• 角度：Θ = 2π U_2 （范围 [0, 2π)）

注意：U_1 必须严格大于 0，否则对数无定义。实际中可采 U_1 ~ U(0,1) 并丢弃 0 值。

3. 计算中间变换的雅可比行列式

第一步：从 (U_1, U_2) 变换到 (R, Θ)。变换关系为：

• R = sqrt(-2 ln U_1) ⇒ U_1 = e^{-R^2/2}
• Θ = 2π U_2 ⇒ U_2 = Θ / (2π)

反向变换的雅可比：计算

$$ J_1 = \left| \frac{\partial (U_1, U_2)}{\partial (R, \Theta)} \right| = \begin{vmatrix} \frac{\partial U_1}{\partial R} & \frac{\partial U_1}{\partial \Theta} \\ \frac{\partial U_2}{\partial R} & \frac{\partial U_2}{\partial \Theta} \end{vmatrix} $$

具体导数：

• ∂U_1/∂R = -R e^{-R^2/2}
• ∂U_1/∂Θ = 0
• ∂U_2/∂R = 0
• ∂U_2/∂Θ = 1/(2π)

因此

$$ J_1 = \left| (-R e^{-R^2/2}) \cdot \frac{1}{2\pi} - 0 \right| = \frac{R}{2\pi} e^{-R^2/2} $$

联合密度：由于 U_1, U_2 独立且密度均为 1（在单位正方形内），联合密度为 f_{U_1,U_2}(u_1,u_2)=1。变换后 (R,Θ) 的联合密度为

$$ f_{R,\Theta}(r,\theta) = 1 \cdot |J_1|^{-1}? $$

需要小心：密度变换公式是 f_{R,Θ}(r,θ) = f_{U_1,U_2}(u_1,u_2) · |J|，其中 J 是从 (U_1,U_2) 到 (R,Θ) 的雅可比，即 ∂(U_1,U_2)/∂(R,Θ)。我们已计算该雅可比为 J_1，所以

$$ f_{R,\Theta}(r,\theta) = 1 \cdot \frac{r}{2\pi} e^{-r^2/2}, \quad r \ge 0, \theta \in [0,2\pi) $$

因子分解：f_{R,Θ} 可写为 (r e^{-r^2/2}) · (1/(2π))，表明 R 和 Θ 独立，且 Θ 服从 U(0,2π)；R 的密度为 f_R(r) = r e^{-r^2/2}（是瑞利分布）。

4. 从 (R,Θ) 变换到笛卡尔坐标 (Z_0, Z_1)

定义：

• Z_0 = R cos Θ
• Z_1 = R sin Θ

逆变换：R = sqrt(Z_0^2 + Z_1^2)，Θ = arctan2(Z_1, Z_0)（需注意象限）。

计算雅可比：从 (R,Θ) 到 (Z_0, Z_1) 的变换是极坐标到直角坐标，雅可比行列式（绝对值）为

$$ J_2 = \left| \frac{\partial (R, \Theta)}{\partial (Z_0, Z_1)} \right| = \frac{1}{R} = \frac{1}{\sqrt{Z_0^2 + Z_1^2}} $$

（因为直角到极坐标的雅可比是 R，逆变换是倒数。）

综合：(Z_0, Z_1) 的联合密度为

$$ f_{Z_0,Z_1}(z_0,z_1) = f_{R,\Theta}(r,\theta) \cdot |J_2|^{-1} ? $$

实际上密度变换公式：f_{Z}(z) = f_{R,Θ}(r(z),θ(z)) · |J_2|，其中 J_2 是从 (R,Θ) 到 (Z_0,Z_1) 的雅可比，但我们已计算其值为 1/R。所以

$$ f_{Z_0,Z_1}(z_0,z_1) = \left(\frac{r}{2\pi} e^{-r^2/2}\right) \cdot \frac{1}{r} = \frac{1}{2\pi} e^{-r^2/2} $$

代入 r^2 = z_0^2 + z_1^2，得

$$ f_{Z_0,Z_1}(z_0,z_1) = \frac{1}{2\pi} e^{-(z_0^2+z_1^2)/2} = \left( \frac{1}{\sqrt{2\pi}} e^{-z_0^2/2} \right) \left( \frac{1}{\sqrt{2\pi}} e^{-z_1^2/2} \right) $$

结论：该联合密度分解为两个标准正态密度乘积，因此 Z_0 和 Z_1 相互独立，且各自服从 N(0,1)。

5. 实现步骤总结

输入：两个独立均匀随机数 u1（0到1，不含0）、u2（0到1）。

执行：

1. 计算 r = sqrt(-2 * ln(u1))
2. 计算 theta = 2 * pi * u2
3. 输出 z0 = r * cos(theta)，z1 = r * sin(theta)

验证：若重复多次，生成的 z0 和 z1 的直方图应近似标准正态曲线，且样本协方差接近 0。

6. 注意数值稳定性和备选形式

• 若 u1 恰好为 0，取对数得到负无穷，实际中应避免。可采 u1 = 1 - U(0,1) 用 -ln(1-u1) 代替。
• 另一种常见形式是“极盒方法”（Marsaglia polar method），避免三角函数调用，思想相同但改用拒绝采样。这里不展开。

最终：现在你已掌握从均匀分布生成正态分布的完整推导。该变换是许多统计软件（如 randn 函数）的底层实现原理。