📚 课程目录第1章均匀分布随机数生成与统计检验

1. 均匀分布随机数生成与统计检验

任务拆解与实现路径

本案例要完成一条完整的“生成—检验—对比”链路：先用线性同余法（LCG）和梅森旋转算法（MT19937）各生成一段 [0,1) 上的均匀随机数序列，再对序列做直方图、K-S 检验和游程检验，最后固定算法、只改变种子，观察序列质量是否稳定。下面给出可直接运行的 Python 实现，并说明每一步的判读依据。

生成器实现

LCG 采用经典参数 a=1664525, c=1013904223, m=2³²，这是 Numerical Recipes 推荐的组合，周期为 2³²。MT19937 直接调用 NumPy 的 default_rng，它内部使用梅森旋转算法，周期远大于 LCG，是工程中的默认选择。

import numpy as np
from scipy import stats

def lcg_uniform(seed, n, a=1664525, c=1013904223, m=2**32):
    """线性同余法生成 [0,1) 均匀随机数"""
    x = seed % m
    out = np.empty(n)
    for i in range(n):
        x = (a * x + c) % m
        out[i] = x / m
    return out

def mt_uniform(seed, n):
    """梅森旋转算法生成 [0,1) 均匀随机数"""
    rng = np.random.default_rng(seed)
    return rng.random(n)

N = 10000
u_lcg = lcg_uniform(seed=42, n=N)
u_mt  = mt_uniform(seed=42, n=N)
print(u_lcg[:5], u_mt[:5])

注意 LCG 的 x 必须用 Python 整数保存，a*x+c 在 2³² 模下不会溢出；若改用 NumPy 的 int32 会溢出出错。MT 版本用 default_rng 而非旧的 np.random.seed，后者是全局状态，不利于种子对比实验。

三项统计检验

直方图：把 [0,1) 等分 10 桶，每桶期望频数 N/10=1000。用卡方拟合优度检验量化偏差，p > 0.05 表示频数分布与均匀分布无显著差异。

K-S 检验：比较经验分布函数 Fₙ(x) 与理论 F(x)=x，统计量 D=max|Fₙ(x)−x|。SciPy 的 kstest(u, 'uniform') 直接给出 D 与 p 值。

游程检验：以序列中位数为阈值，把样本转成 0/1 序列，统计游程数 R。若序列独立，R 的期望与方差为

E[R] = 2*n1*n2/(n1+n2) + 1
Var[R] = 2*n1*n2*(2*n1*n2 - n1 - n2) / ((n1+n2)**2 * (n1+n2-1))
Z = (R - E[R]) / sqrt(Var[R])

|Z| < 1.96 时在 5% 显著性水平下不拒绝独立性假设。

def run_test(u):
    med = np.median(u)
    bits = (u > med).astype(int)
    n1, n2 = bits.sum(), len(bits) - bits.sum()
    runs = 1 + np.sum(bits[1:] != bits[:-1])
    er = 2*n1*n2/(n1+n2) + 1
    var = 2*n1*n2*(2*n1*n2 - n1 - n2) / ((n1+n2)**2 * (n1+n2-1))
    z = (runs - er) / np.sqrt(var)
    return runs, z

def full_check(u, label):
    hist, _ = np.histogram(u, bins=10, range=(0,1))
    chi2, p_chi = stats.chisquare(hist)
    d, p_ks = stats.kstest(u, 'uniform')
    runs, z = run_test(u)
    print(f"{label}: chi2_p={p_chi:.3f}, KS_D={d:.4f}, KS_p={p_ks:.3f}, runs={runs}, Z={z:.3f}")
    return p_chi, p_ks, z

full_check(u_lcg, "LCG seed=42")
full_check(u_mt,  "MT  seed=42")

种子对比实验

固定算法，只改种子，观察三项指标是否稳定。若某个种子下 K-S 的 p 值异常小（如 <0.01）或游程 Z 超出 ±1.96，说明该种子对应的序列段质量偏弱，应换种子或换算法。

for seed in [42, 7, 2024]:
    u = mt_uniform(seed, N)
    full_check(u, f"MT seed={seed}")

for seed in [42, 7, 2024]:
    u = lcg_uniform(seed, N)
    full_check(u, f"LCG seed={seed}")

判读要点：三项检验是“不拒绝”而非“证明均匀”。p 值接近 1 也不代表更好，只要 p > 0.05 且游程 Z 在 ±1.96 内即可认为该序列通过检验。LCG 在低维投影上存在晶格结构，直方图和 K-S 可能通过，但若把相邻两数画成散点图会看到平行线，这是 LCG 的固有缺陷，MT19937 没有这个问题。

常见错误：① 用 np.random.seed 后调用 np.random.rand，种子对比时全局状态被污染；② 游程检验用均值而非中位数分 0/1，导致 n1≠n2 时公式偏差；③ K-S 检验对离散化后的数据使用，会低估 D 值。务必用原始 [0,1) 浮点序列做 K-S。

结果记录表

算法	种子	卡方 p	K-S D	K-S p	游程 Z	结论
LCG	42	—	—	—	—	待运行填写
LCG	7	—	—	—	—	待运行填写
MT19937	42	—	—	—	—	待运行填写
MT19937	2024	—	—	—	—	待运行填写

表中数值需由读者实际运行代码后填入，不同 NumPy/SciPy 版本下 MT19937 的序列一致，LCG 结果完全确定，因此该表可复现。若某行 K-S p < 0.05 或 |Z| > 1.96，标记为“未通过”，并记录该种子，后续案例中避免使用。