4. Box-Muller变换生成二维正态样本

为什么选 Box-Muller：从两个均匀数到一对独立正态数

Box-Muller 变换的核心是：取两个相互独立的均匀随机数 U1, U2 ~ Uniform(0,1)，通过极坐标形式的三角变换，一次性得到两个相互独立的标准正态样本：

Z1 = sqrt(-2·ln U1) · cos(2π·U2)
Z2 = sqrt(-2·ln U1) · sin(2π·U2)

其中 R = sqrt(-2·ln U1) 是半径（服从 Rayleigh 分布），Θ = 2π·U2 是均匀角度。半径与角度独立，因此 Z1 与 Z2 独立且都服从 N(0,1)。这正是它与逆变换法的关键区别：逆变换法一次只产出一个样本，而 Box-Muller 一次产出一对，且天然给出二维联合正态结构，便于直接验证协方差。

实现：一次生成一对样本的向量化写法

下面用 NumPy 实现，注意 U1 必须避开 0（否则 ln 0 = -∞），因此用 np.random.random 后做截断或直接使用 (0,1) 开区间采样。为便于与逆变换法对比，同时给出逆变换法（用 scipy.stats.norm.ppf）的对照实现。

import numpy as np
from scipy import stats

rng = np.random.default_rng(20240501)
N = 200_000                      # 样本对数，最终得到 2N 个正态样本

# ---------- Box-Muller ----------
u1 = rng.random(N)
u2 = rng.random(N)
u1 = np.clip(u1, 1e-12, 1.0)     # 防止 ln(0)
R  = np.sqrt(-2.0 * np.log(u1))
theta = 2.0 * np.pi * u2
z1 = R * np.cos(theta)
z2 = R * np.sin(theta)
bm = np.concatenate([z1, z2])    # 共 2N 个独立标准正态样本

# ---------- 逆变换法（对照） ----------
u = rng.random(2 * N)
inv = stats.norm.ppf(u)          # Φ⁻¹(U)

# ---------- 统计验证 ----------
def report(name, x):
    print(f"{name}: mean={x.mean():+.4f}, var={x.var(ddof=1):.4f}")

report("Box-Muller Z1", z1)
report("Box-Muller Z2", z2)
report("Box-Muller 合并", bm)
report("逆变换法", inv)

# 协方差：Z1 与 Z2 应接近 0
cov = np.cov(z1, z2, ddof=1)[0, 1]
print(f"Cov(Z1,Z2) = {cov:+.5f}")

一次运行（示例种子）的典型输出形态为：均值在 ±0.005 量级、方差在 1.00 附近、协方差在 ±0.005 量级。这些数值随种子变化，属于示例结果，不是固定实测值；判断标准是它们是否落在由样本量决定的合理波动范围内（见下文容差）。

结果判读：均值、方差、协方差与分布对比

对 2N 个样本，均值的标准误约为 1/√(2N)，方差估计的相对标准误约为 √(2/(2N))。以 N=200000 为例，均值容差取 3/√(2N) ≈ 0.0047，方差容差取 3·√(2/(2N)) ≈ 0.0047，协方差容差取 3/√(2N) ≈ 0.0047。若统计量超出该范围，优先检查 U1 是否被截断、角度是否用了 π 而非 2π、以及 cos/sin 是否配对到同一组 R, Θ。

分布对比建议同时画三张 Q-Q 图：Box-Muller 的 Z1、Z2 合并样本，以及逆变换法样本，分别对 N(0,1) 理论分位数作图。Box-Muller 与逆变换法在分布形状上应几乎重合；差异主要体现在计算成本上——Box-Muller 只需 log、sqrt、cos、sin，而逆变换法需要 Φ⁻¹ 的数值近似（如 Acklam 或 scipy 内部实现），单样本成本通常更高。可用 time.perf_counter 对两种方法各跑 10 次取中位数做粗略对比，但需注明这是本机示例计时，不作为通用性能结论。