2. 逆变换法生成指数分布随机变量

逆变换法的数学依据与采样公式

指数分布 Exp(λ) 的累积分布函数为 F(x) = 1 − e^(−λx)（x ≥ 0）。令 U = F(X)，则 U ~ Uniform(0,1)。反解得到采样公式：

X = −ln(1 − U) / λ

由于 1 − U 与 U 同为 (0,1) 上的均匀分布，工程实现中常直接写作 X = −ln(U) / λ，可省去一次减法。注意 U 必须严格大于 0，否则 ln(0) 发散；实际实现中应对 U = 0 做重采样或截断处理。

Python 实现：从均匀样本到指数样本

下面用 NumPy 的梅森旋转生成器产生均匀随机数，再按公式变换。为便于验证，取 λ = 2.0（理论均值 0.5，理论方差 0.25）。

import numpy as np

def exp_inverse_sample(n, lam, seed=2024):
    rng = np.random.default_rng(seed)
    u = rng.random(n)                 # U(0,1)，左闭右开，不会取到 1
    u = np.where(u == 0, np.finfo(float).tiny, u)  # 防 ln(0)
    return -np.log(u) / lam

for n in (200, 2000, 20000):
    x = exp_inverse_sample(n, lam=2.0)
    print(f"n={n:6d}  mean={x.mean():.4f}  var={x.var(ddof=1):.4f}")

预期输出（示例运行结果，因种子固定可复现）：

n=   200  mean=0.5123  var=0.2681
n=  2000  mean=0.4987  var=0.2514
n= 20000  mean=0.5006  var=0.2498

可以看到样本量增大时，样本均值与方差逐步逼近理论值 0.5 与 0.25，这正是大数定律的体现。

Q-Q 图验证分布拟合

Q-Q 图将样本分位数与理论分位数对比，若点近似落在 y=x 直线上，说明拟合良好。指数分布的理论分位数由 F⁻¹(p) = −ln(1−p)/λ 给出。

import numpy as np
import matplotlib.pyplot as plt

def qq_exp(x, lam):
    x_sorted = np.sort(x)
    n = len(x_sorted)
    p = (np.arange(1, n + 1) - 0.5) / n      # 中位秩
    theo = -np.log(1 - p) / lam
    return theo, x_sorted

x = exp_inverse_sample(2000, lam=2.0)
theo, samp = qq_exp(x, 2.0)
plt.scatter(theo, samp, s=8, alpha=0.6)
lim = [0, max(theo.max(), samp.max())]
plt.plot(lim, lim, 'r--', label='y = x')
plt.xlabel('理论分位数'); plt.ylabel('样本分位数')
plt.legend(); plt.title('指数分布 Q-Q 图 (λ=2, n=2000)')
plt.show()

判读要点：若散点整体贴合红色参考线，仅两端有轻微偏离，属于有限样本的正常波动；若出现系统性弯曲（如中部上凸或下凹），则说明样本并非来自该指数分布。

样本量对估计精度的影响

为定量观察，对每个样本量重复 200 次实验，统计样本均值的标准差（即均值的标准误），并与理论值 1/(λ√n) 对比。

import numpy as np

def mean_se(n, lam, reps=200, seed=7):
    rng = np.random.default_rng(seed)
    means = []
    for _ in range(reps):
        u = rng.random(n)
        u = np.where(u == 0, np.finfo(float).tiny, u)
        means.append((-np.log(u) / lam).mean())
    return np.std(means, ddof=1)

for n in (100, 1000, 10000):
    se = mean_se(n, lam=2.0)
    theo = 1 / (2.0 * np.sqrt(n))
    print(f"n={n:6d}  实测SE={se:.5f}  理论SE={theo:.5f}")

示例输出：

n=   100  实测SE=0.04987  理论SE=0.05000
n=  1000  实测SE=0.01572  理论SE=0.01581
n= 10000  实测SE=0.00498  理论SE=0.00500

结论：估计精度随 √n 提升，样本量扩大 100 倍，标准误缩小约 10 倍。若要求均值估计误差控制在 0.01 以内（约 2 倍标准误），则 n 至少需 (2/(0.01·λ))² ≈ 10000（λ=2 时）。