尧图网站建设 
简介压缩包内是一套面向统计信号处理学习者的MATLAB仿真程序围绕广义最大似然比检验GLRT展开可用于研究弱信号检测、噪声背景下异常判别等问题。程序共7个m文件压缩后仅3KB代码体量虽小但模块划分清晰涵盖仿真数据生成、似然比计算、判决门限设定、虚警率与检测概率统计等关键环节适合正在学习假设检验、信号检测理论或需要动手验证算法的本科生与研究生。已有452人学习下载说明该资源对理解似然比检验具有一定参考价值。通过阅读与运行这些脚本读者能直观看到不同信噪比或模型参数下检验性能的变化快速掌握将GLRT理论转化为可运行实验的方法也可基于现有框架扩展自己的检测方案。1. 似然比检验到底是什么从MLE到GLRT一条被低估的统计检验主线做统计检验的人基本都遇到过这种场景手里有一组样本想判断一个参数到底是不是某个值或者有两个嵌套模型想知道复杂模型带来的提升是否显著。用最大似然估计MLE求出参数后接下来最顺理成章的一步就是比较两个模型的最大似然值——这就是似然比检验LRT。当假设里的未知参数被 MLE 替换掉检验就变成广义似然比检验GLRT它是统计检验里覆盖面最广的框架之一。标题里那个spellcw5在复现包里更像是个随机种子标签真正的主角是围绕 GLRT 的构造与实现。这篇笔记会把 GLRT 从理论讲到 Python 落地再讲几个让你翻车的边界坑适合要自己写检验代码、而不是只会调包的从业者。2. 先看懂 GLRT似然比统计量的构造逻辑与适用边界2.1 为什么是“广义”从简单假设到复合假设最原始的似然比检验来自 Neyman-Pearson 引理H0 和 H1 都是简单假设比如 H0: μ0H1: μ1。统计量就是两个似然值的比L(μ0)/L(μ1)比值越大越支持 H0。但真实问题里H1 几乎都是复合假设比如“均值不等于 0”此时 H1 下落在这个集合里的 μ 有无数个你没法直接算一个确定性的似然值。GLRT 的做法很直接把 H0 和 H1 各自参数空间里“最高”的那个点拿出来比较。也就是先在 H0 约束下找最大似然估计再在无约束全集参数空间里找最大似然估计然后看这两个最高峰的比值。因为没有直接指定参数值而是用 MLE 去估计所以叫“广义”似然比检验。这个“广义”两个字恰恰是它好用的原因——把复合假设压缩成两个最高点的比较问题一下子就可算了。形式化写法是H0 下的最大对数似然l0 sup_{θ∈Θ0} log L(θ)无约束的最大对数似然l1 sup_{θ∈Θ} log L(θ)检验统计量LR 2 * (l1 - l0)为什么前面乘 2这是为了让它渐近服从卡方分布。Wilks 定理在 1938 年给出了结论在 H0 成立且满足一系列正则条件时LR会依分布收敛到 χ²自由度等于约束掉参数的个数。这个结论是整个 GLRT 能落地的基石也是后面很多坑的源头。2.2 MLE 在 GLRT 中的角色谁来做无约束和约束估计很多人用 GLRT 时会忽略一个关键点MLE 不是被检验的对象而是构造统计量的工具。无约束 MLE 在全集参数空间找最高峰约束 MLE 在 H0 定义的低维子空间里找最高峰。两个峰的高度差就是数据对 H0 的“意见”。如果约束 MLE 恰好等于无约束 MLE比如 H0 里要求的均值正是样本均值那么 LR0p 值会是 1。这种情况不是“数据强烈支持 H0”而是“数据连区分 H0 和 H1 的能力都没有”。你只能说在这个样本下两个模型没有差异能被检测出来。MLE 的求解有两种路径有闭式解的正态分布直接用公式没有闭式解的模型比如逻辑回归、混合模型就得用数值优化。数值优化时要注意约束 MLE 并不是简单地“把无约束 MLE 的某些参数设成固定值”再重新优化而可能需要在约束流形上重新搜索。有些约束会让参数空间边界出现这就是后面第 4 章的坑。2.3 渐近分布与自由度卡方近似的使用条件教科书里通常写LR ~ χ²(q)其中 q 是自由度等于无约束参数个数减去约束参数个数。这句话有个前置条件真参数必须在参数空间内部模型可识别且样本量足够大。这三个条件任何一个被破坏卡方近似就会失效。举两个典型的失效场景。第一个是边界参数比如检验 H0: σ² ≤ σ0²如果约束解恰好落在边界 σ0² 上那么 LR 的渐近分布不再是 χ²(1)而是一个 50:50 的混合卡方。第二个是模型不可识别混合模型里检验稀疏成分是否存在时H0 下某些参数失去辨识度LR 的渐近分布可能是非标准分布。所以使用卡方近似前先看一眼参数空间的结构会比直接算 p 值安全得多。3. 用 Python 实现 GLRT完整复现一个正态均值检验3.1 最小可跑通代码生成数据并计算 LR 统计量我一般会用最简单的正态分布均值检验来验证我对 GLRT 的理解有没有跑偏。假设数据来自N(μ, σ²)方差也未知检验目标是 H0: μ0 对 H1: μ≠0。这个模型有两个参数μ, σ²约束条件是 μ0所以自由度是 1。正态分布的 MLE 有闭式解不需要开优化器无约束μ_hat mean(x)σ²_hat mean((x - μ_hat)^2)约束μ0σ²_0 mean(x^2)两个模型的最大对数似然之差可以化简成LR n * log(σ²_0 / σ²_hat)。代码实现如下import numpy as np from scipy import stats # 固定随机种子确保结果可复现 rng np.random.default_rng(20240617) n 60 mu_true 0.4 sigma_true 1.0 x rng.normal(mu_true, sigma_true, sizen) def fit_gaussian(data, mu_fixedNone): 返回高斯MLE的(mu, sigma2) 当mu_fixed给定时只估计sigma2约束模型 if mu_fixed is None: mu data.mean() else: mu mu_fixed sigma2 np.mean((data - mu) ** 2) return mu, sigma2 # 无约束MLE mu_mle, sigma2_mle fit_gaussian(x) # 约束MLEH0: mu0 mu0, sigma2_0 fit_gaussian(x, mu_fixed0.0) # 对数似然比统计量 lr n * np.log(sigma2_0 / sigma2_mle) df 1 p_value stats.chi2.sf(lr, df) print(f无约束MLE: mu{mu_mle:.4f}, sigma2{sigma2_mle:.4f}) print(f约束MLE: mu0, sigma2{sigma2_0:.4f}) print(fLR statistic: {lr:.4f}) print(fp-value (chi2 df1): {p_value:.4f})代码里的fit_gaussian返回的是以 n 为分母的方差估计这才是最大似然估计如果用 n-1 做分母算出来的广义似然比会和真正定义不一致。stats.chi2.sf(lr, df)算的是上尾概率也就是 p 值这比1 - stats.chi2.cdf(lr, df)更稳能避免浮点数在尾部产生负零。这里df1是因为 H0 只约束了 μ 一个参数σ² 一直在被估计。3.2 关键参数怎么设样本量、方差估计与 p 值计算实际落地时我最常被问到的就是 n 到底取多少。GLRT 是渐近检验n 越大卡方近似越准但“大到多少”没有绝对标准。经验上n 小于 30 时建议用 Bootstrap 做一次校准再下结论n 在 50 以上如果参数不在边界χ² 近似基本可信。另外效应量越小需要的样本量越大。如果真实效应量只有 0.1 个标准差n30 的 GLRT 功效会低得可怜p 值不显著不代表没有效应。方差估计用 n 还是 n-1这里不是随意选的。因为对数似然定义里的方差就是mean((x-μ)^2)如果换成无偏估计var(x)LR 统计量就不是严格意义上的似然比。虽然数值差距通常很小但在写论文或交付代码时定义一定要和公式对上。显著性水平 α 默认 0.05但如果做多重检验需要先做 Bonferroni 或 FDR 校正再套用chi2.sf的临界值。不要拿校正后的 p 值反推 LR 阈值那是另一套逻辑。3.3 从闭式解到通用优化器当模型不再有解析 MLE正态模型的闭式解帮你验证了原理但大多数模型没有这种好运。逻辑回归、威布尔分布、混合模型都得靠数值优化。下面这段代码展示怎么用scipy.optimize.minimize做无约束和约束最大化。from scipy.optimize import minimize def neg_loglik_gaussian(theta, data): mu, log_sigma2 theta sigma2 np.exp(log_sigma2) n len(data) # 对数似然去掉与mu、sigma2无关的常数项 ll -0.5 * np.sum((data - mu) ** 2) / sigma2 - 0.5 * n * np.log(sigma2) return -ll # 无约束优化 theta0 np.array([0.0, 0.0]) res1 minimize(neg_loglik_gaussian, theta0, args(x,), methodL-BFGS-B) mu_mle res1.x[0] sigma2_mle np.exp(res1.x[1]) l1 -res1.fun # 约束优化固定mu0只优化log_sigma2 def neg_loglik_constrained(log_sigma2, data): return neg_loglik_gaussian(np.array([0.0, log_sigma2[0]]), data) res0 minimize(neg_loglik_constrained, np.array([0.0]), args(x,), methodL-BFGS-B) sigma2_0 np.exp(res0.x[0]) l0 -res0.fun lr_opt 2 * (l1 - l0) p_opt stats.chi2.sf(lr_opt, 1) print(f优化版 LR: {lr_opt:.4f}, p: {p_opt:.4f})注意这段代码里我把方差参数化成了log_sigma2。原因很简单方差必须为正直接优化 σ² 很容易让优化器闯进负数区域导致似然函数无定义。用 log 变换后参数空间变成整个实数轴无约束优化就不再受边界限制。L-BFGS-B适合中小规模光滑优化问题比BFGS更抗噪声。如果你发现优化结果和闭式解差很远先检查是不是初始值给得太离谱我通常会从矩估计出发而不是拍脑袋给 0。4. GLRT 常见问题与排查4 个典型坑从“拒绝”到“错误拒绝”的距离4.1 边界参数让卡方近似失效现象你做了 GLRT算出来的 p 值明显和模拟不符。比如 H0: σ² ≤ σ0²样本方差刚好略低于 σ0²LR 很小但用 χ²(1) 算出的 p 值总是偏大。原因当约束估计落在参数空间边界上时Wilks 定理的正则条件不成立。H0 的真值在边界Fisher 信息矩阵的渐近特性变了LR 的渐近分布是一个 50:50 的混合一半概率为 0一半概率为 χ²(1)。解决这类检验不能直接套chi2.sf。正确的 p 值应该用0.5 0.5 * stats.chi2.sf(lr, 1)或更一般的混合分布计算。如果你不确定边界情况最稳妥的办法是跑一次 Bootstrap让数据自己告诉你分布长什么样。我过去吃过一次亏用标准 χ² 给边界检验算 p 值结果在一个仿真任务里拒绝了 H0后来换 Bootstrap 才发现根本不显著。4.2 似然函数写错一个符号统计量变成负的现象LR 算出来是负数p 值大于 1。这种情况在数值优化里很常见尤其是手工推导对数似然时漏掉了某个常数项。原因无约束最大对数似然理论上一定大于等于约束最大对数似然。如果出现l1 l0几乎必然是似然函数写错。最常见的是少了-0.5 * n * log(sigma2)这一项或者把平方和的负号漏掉。解决先用最简单的闭式解做一个 sanity check。比如在第 3.1 节的代码里打印l1和l0如果l1 l0立刻检查对数似然公式。另外优化器返回的res.fun是负对数似然记得取相反数。很多初学踩坑的人就是忘了l1 -res1.fun直接把res1.fun拿来当正值。4.3 约束估计不收敛换算法还是换初始值现象minimize返回successFalse或者警告Desired error not necessarily achieved due to precision loss。迭代次数打满结果却没落在约束附近。原因统计优化里九成的不收敛都出在初始值。目标函数在高维空间可能存在平坦区域或者不同参数尺度差异太大比如 μ 在 0.1 量级而 σ² 在 1000 量级优化器会来回震荡。解决先换初始值。不要从零开始用矩估计作为起点mu0 x.mean()sigma2_0 x.var()。如果还不行就做参数重参数化用 log 变换把尺度拉齐。最后才是换优化算法Nelder-Mead对初始值不敏感但收敛慢L-BFGS-B快但对光滑性有要求。另外约束优化不要自己硬编码惩罚项直接用scipy.optimize的constraints参数或固定参数后降维比我早期手动加penalty可靠得多。4.4 小样本下 p 值偏保守先看功效曲线再下结论现象n20 时χ² 近似给出的 p 值是 0.08但参数 Bootstrap 给出的 p 值是 0.03。一个不显著一个显著。原因小样本下 LR 统计量的分布并不是 χ²而是更集中在左侧也就是说尾部比 χ² 更薄。用 χ² 的上尾概率会高估 p 值导致检验偏保守。保守听起来不危险但在样本量受限的医学或工业检测里它会让真实的效应被当成噪声放过。解决小样本下优先用 Bootstrap 校准至少作为敏感性分析上报。另一个习惯是画出功效曲线看看在当前样本量下这个检验对多大的效应量有 80% 的把握。如果效应量本身就低于能检出的阈值p 值偏大并不意外。保守不是 bug但要清楚它有多保守。4.5 模型不可识别的时候自由度不是参数个数差现象检验混合模型里是否存在两个成分H0 设定两个成分均值相等。你用参数个数差当自由度算出 LR 对应一个很宽的 p 值结果和仿真对不上。原因当 H0 成立时混合比例这个参数根本不可辨识——两个均值相同的成分等价于一个成分混合比例可以被任意值替代。Fisher 信息矩阵退化Wilks 定理不适用LR 的渐近分布不是 χ²。解决不要手写自由度。要么做参数 Bootstrap要么查专门针对不可识别模型的检验理论比如 Chernoff 给出的非标准极限分布。实操中Bootstrap 最省事。我在隐马尔可夫模型的项目里遇到过同样的问题最后是换成仿真 p 值才解决了评审质疑。这类问题没有通用自由度公式数据仿真永远比拍脑袋安全。5. 理论分布靠不住时用 Bootstrap 校准 GLRT 的 p 值5.1 什么时候必须放弃卡方近似先做一次 Bootstrap 对比判断标准不是“n 小于多少”而是看看 χ² 近似在你这个模型下到底靠不靠谱。我有一个低成本的做法先跑 200 次参数 Bootstrap把 LR 经验分布和 χ² 曲线叠在一起看。如果两条线在尾部明显分离就说明近似失效。触发分离的常见条件有三个一是样本量小尤其是 n 小于 30二是参数落在边界比如约束方差、约束概率为 0三是模型存在不可识别参数比如混合模型、随机效应方差为 0 的情况。这些条件下χ² 的尾部可能偏厚也可能偏薄方向不确定所以不要猜去仿。5.2 参数 Bootstrap 实现在 H0 下生成数据重新计算 LR参数 Bootstrap 的核心思想是从拟合好的 H0 模型里生成样本再对每个样本重新做一遍 GLRT从而得到 H0 下 LR 的经验分布。下面的代码接第 3.1 节的正态检验例子。B 999 # 基于原始数据在 H0 下的约束估计 _, sigma2_0 fit_gaussian(x, mu_fixed0.0) lr_boot np.empty(B) for b in range(B): # 从 H0 模型生成数据 xb rng.normal(0.0, np.sqrt(sigma2_0), sizen) # 对每个 bootstrap 样本重新估计两个模型 mu_mle_b xb.mean() sigma2_mle_b np.mean((xb - mu_mle_b) ** 2) sigma2_0_b np.mean(xb ** 2) lr_b n * np.log(sigma2_0_b / sigma2_mle_b) lr_boot[b] lr_b p_boot (1 np.sum(lr_boot lr)) / (B 1) print(fBootstrap p-value: {p_boot:.4f})生成数据时用的是原始数据在 H0 下的约束估计sigma2_0保证样本确实来自 H0。注意lr是原始数据的观测统计量Bootstrap 循环里重新算的是每个仿真样本自己的统计量。p 值公式里分子加 1、分母加 1是为了防止 p 值变成 0也保证有限样本下的有效 p 值不会低于 1/(B1)。如果你要更高的精度B 可以取 9999但要注意时间成本。这段代码里有一个隐藏细节每个 Bootstrap 样本的约束估计sigma2_0_b是重新从xb算的而不是沿用原始数据的sigma2_0。否则你会把原始样本的信息带进 H0 分布Bootstrap 就变成半参数版本了结果会偏乐观。5.3 怎么解读校准结果p 值变保守还是激进拿到 Bootstrap p 值和 χ² p 值后对比一下方向。如果 Bootstrap p 更大说明 χ² 近似在往激进的方向偏也就是更容易拒绝 H0这时再用 χ² 会有假阳性风险。如果 Bootstrap p 更小说明 χ² 保守可能漏检。不要只说“Bootstrap 更准”要报告两者差异有多大。我见过不少人在论文里只写 Bootstrap p 值不解释为什么不用 χ²。稳妥的做法是两种 p 值都报告并注明“由于样本量较小/参数位于边界采用参数 Bootstrap 校准”。这种透明表述比任何检验技巧都能增加可信度。另外Bootstrap 结果受随机数种子影响务必固定种子否则别人复现你的 p 值时会对不上。6. 一个能救命的习惯把 GLRT 的拒绝域画出来只报 p 值远远不够。我最近两年养成了一个习惯每次用 GLRT 之前先模拟出它的拒绝域或功效曲线。这一步能直接暴露样本量和效应量的问题比事后盯着 p 值猜原因有用得多。下面这段代码模拟了不同样本量和不同真实均值下GLRT 的拒绝概率。注意这里用了第 3.1 节的方差未知 LR 公式但人为固定了方差为 1 来简化模拟。def glrt_reject(mu_true, n, alpha0.05, R500): rej 0 for _ in range(R): xsim rng.normal(mu_true, 1.0, n) lr_sim n * np.log(np.mean(xsim**2) / np.mean((xsim - xsim.mean())**2)) if stats.chi2.sf(lr_sim, 1) alpha: rej 1 return rej / R mu_grid np.linspace(0, 1.2, 13) power_n20 [glrt_reject(m, 20) for m in mu_grid] power_n100 [glrt_reject(m, 100) for m in mu_grid] # 画图保持简洁略去图形样式 import matplotlib.pyplot as plt plt.plot(mu_grid, power_n20, markero, labeln20) plt.plot(mu_grid, power_n100, markers, labeln100) plt.xlabel(true mean mu) plt.ylabel(probability of rejecting H0) plt.legend() plt.show()从图上你会看到一条很扎眼的规律真实效应量为 0.4 时n20 的检验可能只有 40% 左右的把握拒绝 H0而 n100 能达到 90% 以上。也就是说如果你在 n20 时看到 p0.12先别急着写“没有效应”很可能是检验根本没带耳朵。这个模拟图就是你的后悔药在投入数据采集之前先问自己“我要找的效应量在这个样本量下检不检得出”。我的一个血泪经验是有一年用 GLRT 做特征选择只看了 p 值结果重要特征因为样本量不足被过滤掉。后来画出功效曲线才发现那个效应量在当时的 n 下只有 30% 的把握能检出来。从那以后我给自己定了条规矩任何 GLRT 结论必须附带一组功效曲线或者至少说明最小可检测效应量。这个习惯看起来多花十分钟但能避免你对着一个没有分辨能力的统计量瞎解释。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取 GLRT似然比检验原理与Python实现：从MLE到卡方近似及Bootstrap校准
简介压缩包内是一套面向统计信号处理学习者的MATLAB仿真程序围绕广义最大似然比检验GLRT展开可用于研究弱信号检测、噪声背景下异常判别等问题。程序共7个m文件压缩后仅3KB代码体量虽小但模块划分清晰涵盖仿真数据生成、似然比计算、判决门限设定、虚警率与检测概率统计等关键环节适合正在学习假设检验、信号检测理论或需要动手验证算法的本科生与研究生。已有452人学习下载说明该资源对理解似然比检验具有一定参考价值。通过阅读与运行这些脚本读者能直观看到不同信噪比或模型参数下检验性能的变化快速掌握将GLRT理论转化为可运行实验的方法也可基于现有框架扩展自己的检测方案。1. 似然比检验到底是什么从MLE到GLRT一条被低估的统计检验主线做统计检验的人基本都遇到过这种场景手里有一组样本想判断一个参数到底是不是某个值或者有两个嵌套模型想知道复杂模型带来的提升是否显著。用最大似然估计MLE求出参数后接下来最顺理成章的一步就是比较两个模型的最大似然值——这就是似然比检验LRT。当假设里的未知参数被 MLE 替换掉检验就变成广义似然比检验GLRT它是统计检验里覆盖面最广的框架之一。标题里那个spellcw5在复现包里更像是个随机种子标签真正的主角是围绕 GLRT 的构造与实现。这篇笔记会把 GLRT 从理论讲到 Python 落地再讲几个让你翻车的边界坑适合要自己写检验代码、而不是只会调包的从业者。2. 先看懂 GLRT似然比统计量的构造逻辑与适用边界2.1 为什么是“广义”从简单假设到复合假设最原始的似然比检验来自 Neyman-Pearson 引理H0 和 H1 都是简单假设比如 H0: μ0H1: μ1。统计量就是两个似然值的比L(μ0)/L(μ1)比值越大越支持 H0。但真实问题里H1 几乎都是复合假设比如“均值不等于 0”此时 H1 下落在这个集合里的 μ 有无数个你没法直接算一个确定性的似然值。GLRT 的做法很直接把 H0 和 H1 各自参数空间里“最高”的那个点拿出来比较。也就是先在 H0 约束下找最大似然估计再在无约束全集参数空间里找最大似然估计然后看这两个最高峰的比值。因为没有直接指定参数值而是用 MLE 去估计所以叫“广义”似然比检验。这个“广义”两个字恰恰是它好用的原因——把复合假设压缩成两个最高点的比较问题一下子就可算了。形式化写法是H0 下的最大对数似然l0 sup_{θ∈Θ0} log L(θ)无约束的最大对数似然l1 sup_{θ∈Θ} log L(θ)检验统计量LR 2 * (l1 - l0)为什么前面乘 2这是为了让它渐近服从卡方分布。Wilks 定理在 1938 年给出了结论在 H0 成立且满足一系列正则条件时LR会依分布收敛到 χ²自由度等于约束掉参数的个数。这个结论是整个 GLRT 能落地的基石也是后面很多坑的源头。2.2 MLE 在 GLRT 中的角色谁来做无约束和约束估计很多人用 GLRT 时会忽略一个关键点MLE 不是被检验的对象而是构造统计量的工具。无约束 MLE 在全集参数空间找最高峰约束 MLE 在 H0 定义的低维子空间里找最高峰。两个峰的高度差就是数据对 H0 的“意见”。如果约束 MLE 恰好等于无约束 MLE比如 H0 里要求的均值正是样本均值那么 LR0p 值会是 1。这种情况不是“数据强烈支持 H0”而是“数据连区分 H0 和 H1 的能力都没有”。你只能说在这个样本下两个模型没有差异能被检测出来。MLE 的求解有两种路径有闭式解的正态分布直接用公式没有闭式解的模型比如逻辑回归、混合模型就得用数值优化。数值优化时要注意约束 MLE 并不是简单地“把无约束 MLE 的某些参数设成固定值”再重新优化而可能需要在约束流形上重新搜索。有些约束会让参数空间边界出现这就是后面第 4 章的坑。2.3 渐近分布与自由度卡方近似的使用条件教科书里通常写LR ~ χ²(q)其中 q 是自由度等于无约束参数个数减去约束参数个数。这句话有个前置条件真参数必须在参数空间内部模型可识别且样本量足够大。这三个条件任何一个被破坏卡方近似就会失效。举两个典型的失效场景。第一个是边界参数比如检验 H0: σ² ≤ σ0²如果约束解恰好落在边界 σ0² 上那么 LR 的渐近分布不再是 χ²(1)而是一个 50:50 的混合卡方。第二个是模型不可识别混合模型里检验稀疏成分是否存在时H0 下某些参数失去辨识度LR 的渐近分布可能是非标准分布。所以使用卡方近似前先看一眼参数空间的结构会比直接算 p 值安全得多。3. 用 Python 实现 GLRT完整复现一个正态均值检验3.1 最小可跑通代码生成数据并计算 LR 统计量我一般会用最简单的正态分布均值检验来验证我对 GLRT 的理解有没有跑偏。假设数据来自N(μ, σ²)方差也未知检验目标是 H0: μ0 对 H1: μ≠0。这个模型有两个参数μ, σ²约束条件是 μ0所以自由度是 1。正态分布的 MLE 有闭式解不需要开优化器无约束μ_hat mean(x)σ²_hat mean((x - μ_hat)^2)约束μ0σ²_0 mean(x^2)两个模型的最大对数似然之差可以化简成LR n * log(σ²_0 / σ²_hat)。代码实现如下import numpy as np from scipy import stats # 固定随机种子确保结果可复现 rng np.random.default_rng(20240617) n 60 mu_true 0.4 sigma_true 1.0 x rng.normal(mu_true, sigma_true, sizen) def fit_gaussian(data, mu_fixedNone): 返回高斯MLE的(mu, sigma2) 当mu_fixed给定时只估计sigma2约束模型 if mu_fixed is None: mu data.mean() else: mu mu_fixed sigma2 np.mean((data - mu) ** 2) return mu, sigma2 # 无约束MLE mu_mle, sigma2_mle fit_gaussian(x) # 约束MLEH0: mu0 mu0, sigma2_0 fit_gaussian(x, mu_fixed0.0) # 对数似然比统计量 lr n * np.log(sigma2_0 / sigma2_mle) df 1 p_value stats.chi2.sf(lr, df) print(f无约束MLE: mu{mu_mle:.4f}, sigma2{sigma2_mle:.4f}) print(f约束MLE: mu0, sigma2{sigma2_0:.4f}) print(fLR statistic: {lr:.4f}) print(fp-value (chi2 df1): {p_value:.4f})代码里的fit_gaussian返回的是以 n 为分母的方差估计这才是最大似然估计如果用 n-1 做分母算出来的广义似然比会和真正定义不一致。stats.chi2.sf(lr, df)算的是上尾概率也就是 p 值这比1 - stats.chi2.cdf(lr, df)更稳能避免浮点数在尾部产生负零。这里df1是因为 H0 只约束了 μ 一个参数σ² 一直在被估计。3.2 关键参数怎么设样本量、方差估计与 p 值计算实际落地时我最常被问到的就是 n 到底取多少。GLRT 是渐近检验n 越大卡方近似越准但“大到多少”没有绝对标准。经验上n 小于 30 时建议用 Bootstrap 做一次校准再下结论n 在 50 以上如果参数不在边界χ² 近似基本可信。另外效应量越小需要的样本量越大。如果真实效应量只有 0.1 个标准差n30 的 GLRT 功效会低得可怜p 值不显著不代表没有效应。方差估计用 n 还是 n-1这里不是随意选的。因为对数似然定义里的方差就是mean((x-μ)^2)如果换成无偏估计var(x)LR 统计量就不是严格意义上的似然比。虽然数值差距通常很小但在写论文或交付代码时定义一定要和公式对上。显著性水平 α 默认 0.05但如果做多重检验需要先做 Bonferroni 或 FDR 校正再套用chi2.sf的临界值。不要拿校正后的 p 值反推 LR 阈值那是另一套逻辑。3.3 从闭式解到通用优化器当模型不再有解析 MLE正态模型的闭式解帮你验证了原理但大多数模型没有这种好运。逻辑回归、威布尔分布、混合模型都得靠数值优化。下面这段代码展示怎么用scipy.optimize.minimize做无约束和约束最大化。from scipy.optimize import minimize def neg_loglik_gaussian(theta, data): mu, log_sigma2 theta sigma2 np.exp(log_sigma2) n len(data) # 对数似然去掉与mu、sigma2无关的常数项 ll -0.5 * np.sum((data - mu) ** 2) / sigma2 - 0.5 * n * np.log(sigma2) return -ll # 无约束优化 theta0 np.array([0.0, 0.0]) res1 minimize(neg_loglik_gaussian, theta0, args(x,), methodL-BFGS-B) mu_mle res1.x[0] sigma2_mle np.exp(res1.x[1]) l1 -res1.fun # 约束优化固定mu0只优化log_sigma2 def neg_loglik_constrained(log_sigma2, data): return neg_loglik_gaussian(np.array([0.0, log_sigma2[0]]), data) res0 minimize(neg_loglik_constrained, np.array([0.0]), args(x,), methodL-BFGS-B) sigma2_0 np.exp(res0.x[0]) l0 -res0.fun lr_opt 2 * (l1 - l0) p_opt stats.chi2.sf(lr_opt, 1) print(f优化版 LR: {lr_opt:.4f}, p: {p_opt:.4f})注意这段代码里我把方差参数化成了log_sigma2。原因很简单方差必须为正直接优化 σ² 很容易让优化器闯进负数区域导致似然函数无定义。用 log 变换后参数空间变成整个实数轴无约束优化就不再受边界限制。L-BFGS-B适合中小规模光滑优化问题比BFGS更抗噪声。如果你发现优化结果和闭式解差很远先检查是不是初始值给得太离谱我通常会从矩估计出发而不是拍脑袋给 0。4. GLRT 常见问题与排查4 个典型坑从“拒绝”到“错误拒绝”的距离4.1 边界参数让卡方近似失效现象你做了 GLRT算出来的 p 值明显和模拟不符。比如 H0: σ² ≤ σ0²样本方差刚好略低于 σ0²LR 很小但用 χ²(1) 算出的 p 值总是偏大。原因当约束估计落在参数空间边界上时Wilks 定理的正则条件不成立。H0 的真值在边界Fisher 信息矩阵的渐近特性变了LR 的渐近分布是一个 50:50 的混合一半概率为 0一半概率为 χ²(1)。解决这类检验不能直接套chi2.sf。正确的 p 值应该用0.5 0.5 * stats.chi2.sf(lr, 1)或更一般的混合分布计算。如果你不确定边界情况最稳妥的办法是跑一次 Bootstrap让数据自己告诉你分布长什么样。我过去吃过一次亏用标准 χ² 给边界检验算 p 值结果在一个仿真任务里拒绝了 H0后来换 Bootstrap 才发现根本不显著。4.2 似然函数写错一个符号统计量变成负的现象LR 算出来是负数p 值大于 1。这种情况在数值优化里很常见尤其是手工推导对数似然时漏掉了某个常数项。原因无约束最大对数似然理论上一定大于等于约束最大对数似然。如果出现l1 l0几乎必然是似然函数写错。最常见的是少了-0.5 * n * log(sigma2)这一项或者把平方和的负号漏掉。解决先用最简单的闭式解做一个 sanity check。比如在第 3.1 节的代码里打印l1和l0如果l1 l0立刻检查对数似然公式。另外优化器返回的res.fun是负对数似然记得取相反数。很多初学踩坑的人就是忘了l1 -res1.fun直接把res1.fun拿来当正值。4.3 约束估计不收敛换算法还是换初始值现象minimize返回successFalse或者警告Desired error not necessarily achieved due to precision loss。迭代次数打满结果却没落在约束附近。原因统计优化里九成的不收敛都出在初始值。目标函数在高维空间可能存在平坦区域或者不同参数尺度差异太大比如 μ 在 0.1 量级而 σ² 在 1000 量级优化器会来回震荡。解决先换初始值。不要从零开始用矩估计作为起点mu0 x.mean()sigma2_0 x.var()。如果还不行就做参数重参数化用 log 变换把尺度拉齐。最后才是换优化算法Nelder-Mead对初始值不敏感但收敛慢L-BFGS-B快但对光滑性有要求。另外约束优化不要自己硬编码惩罚项直接用scipy.optimize的constraints参数或固定参数后降维比我早期手动加penalty可靠得多。4.4 小样本下 p 值偏保守先看功效曲线再下结论现象n20 时χ² 近似给出的 p 值是 0.08但参数 Bootstrap 给出的 p 值是 0.03。一个不显著一个显著。原因小样本下 LR 统计量的分布并不是 χ²而是更集中在左侧也就是说尾部比 χ² 更薄。用 χ² 的上尾概率会高估 p 值导致检验偏保守。保守听起来不危险但在样本量受限的医学或工业检测里它会让真实的效应被当成噪声放过。解决小样本下优先用 Bootstrap 校准至少作为敏感性分析上报。另一个习惯是画出功效曲线看看在当前样本量下这个检验对多大的效应量有 80% 的把握。如果效应量本身就低于能检出的阈值p 值偏大并不意外。保守不是 bug但要清楚它有多保守。4.5 模型不可识别的时候自由度不是参数个数差现象检验混合模型里是否存在两个成分H0 设定两个成分均值相等。你用参数个数差当自由度算出 LR 对应一个很宽的 p 值结果和仿真对不上。原因当 H0 成立时混合比例这个参数根本不可辨识——两个均值相同的成分等价于一个成分混合比例可以被任意值替代。Fisher 信息矩阵退化Wilks 定理不适用LR 的渐近分布不是 χ²。解决不要手写自由度。要么做参数 Bootstrap要么查专门针对不可识别模型的检验理论比如 Chernoff 给出的非标准极限分布。实操中Bootstrap 最省事。我在隐马尔可夫模型的项目里遇到过同样的问题最后是换成仿真 p 值才解决了评审质疑。这类问题没有通用自由度公式数据仿真永远比拍脑袋安全。5. 理论分布靠不住时用 Bootstrap 校准 GLRT 的 p 值5.1 什么时候必须放弃卡方近似先做一次 Bootstrap 对比判断标准不是“n 小于多少”而是看看 χ² 近似在你这个模型下到底靠不靠谱。我有一个低成本的做法先跑 200 次参数 Bootstrap把 LR 经验分布和 χ² 曲线叠在一起看。如果两条线在尾部明显分离就说明近似失效。触发分离的常见条件有三个一是样本量小尤其是 n 小于 30二是参数落在边界比如约束方差、约束概率为 0三是模型存在不可识别参数比如混合模型、随机效应方差为 0 的情况。这些条件下χ² 的尾部可能偏厚也可能偏薄方向不确定所以不要猜去仿。5.2 参数 Bootstrap 实现在 H0 下生成数据重新计算 LR参数 Bootstrap 的核心思想是从拟合好的 H0 模型里生成样本再对每个样本重新做一遍 GLRT从而得到 H0 下 LR 的经验分布。下面的代码接第 3.1 节的正态检验例子。B 999 # 基于原始数据在 H0 下的约束估计 _, sigma2_0 fit_gaussian(x, mu_fixed0.0) lr_boot np.empty(B) for b in range(B): # 从 H0 模型生成数据 xb rng.normal(0.0, np.sqrt(sigma2_0), sizen) # 对每个 bootstrap 样本重新估计两个模型 mu_mle_b xb.mean() sigma2_mle_b np.mean((xb - mu_mle_b) ** 2) sigma2_0_b np.mean(xb ** 2) lr_b n * np.log(sigma2_0_b / sigma2_mle_b) lr_boot[b] lr_b p_boot (1 np.sum(lr_boot lr)) / (B 1) print(fBootstrap p-value: {p_boot:.4f})生成数据时用的是原始数据在 H0 下的约束估计sigma2_0保证样本确实来自 H0。注意lr是原始数据的观测统计量Bootstrap 循环里重新算的是每个仿真样本自己的统计量。p 值公式里分子加 1、分母加 1是为了防止 p 值变成 0也保证有限样本下的有效 p 值不会低于 1/(B1)。如果你要更高的精度B 可以取 9999但要注意时间成本。这段代码里有一个隐藏细节每个 Bootstrap 样本的约束估计sigma2_0_b是重新从xb算的而不是沿用原始数据的sigma2_0。否则你会把原始样本的信息带进 H0 分布Bootstrap 就变成半参数版本了结果会偏乐观。5.3 怎么解读校准结果p 值变保守还是激进拿到 Bootstrap p 值和 χ² p 值后对比一下方向。如果 Bootstrap p 更大说明 χ² 近似在往激进的方向偏也就是更容易拒绝 H0这时再用 χ² 会有假阳性风险。如果 Bootstrap p 更小说明 χ² 保守可能漏检。不要只说“Bootstrap 更准”要报告两者差异有多大。我见过不少人在论文里只写 Bootstrap p 值不解释为什么不用 χ²。稳妥的做法是两种 p 值都报告并注明“由于样本量较小/参数位于边界采用参数 Bootstrap 校准”。这种透明表述比任何检验技巧都能增加可信度。另外Bootstrap 结果受随机数种子影响务必固定种子否则别人复现你的 p 值时会对不上。6. 一个能救命的习惯把 GLRT 的拒绝域画出来只报 p 值远远不够。我最近两年养成了一个习惯每次用 GLRT 之前先模拟出它的拒绝域或功效曲线。这一步能直接暴露样本量和效应量的问题比事后盯着 p 值猜原因有用得多。下面这段代码模拟了不同样本量和不同真实均值下GLRT 的拒绝概率。注意这里用了第 3.1 节的方差未知 LR 公式但人为固定了方差为 1 来简化模拟。def glrt_reject(mu_true, n, alpha0.05, R500): rej 0 for _ in range(R): xsim rng.normal(mu_true, 1.0, n) lr_sim n * np.log(np.mean(xsim**2) / np.mean((xsim - xsim.mean())**2)) if stats.chi2.sf(lr_sim, 1) alpha: rej 1 return rej / R mu_grid np.linspace(0, 1.2, 13) power_n20 [glrt_reject(m, 20) for m in mu_grid] power_n100 [glrt_reject(m, 100) for m in mu_grid] # 画图保持简洁略去图形样式 import matplotlib.pyplot as plt plt.plot(mu_grid, power_n20, markero, labeln20) plt.plot(mu_grid, power_n100, markers, labeln100) plt.xlabel(true mean mu) plt.ylabel(probability of rejecting H0) plt.legend() plt.show()从图上你会看到一条很扎眼的规律真实效应量为 0.4 时n20 的检验可能只有 40% 左右的把握拒绝 H0而 n100 能达到 90% 以上。也就是说如果你在 n20 时看到 p0.12先别急着写“没有效应”很可能是检验根本没带耳朵。这个模拟图就是你的后悔药在投入数据采集之前先问自己“我要找的效应量在这个样本量下检不检得出”。我的一个血泪经验是有一年用 GLRT 做特征选择只看了 p 值结果重要特征因为样本量不足被过滤掉。后来画出功效曲线才发现那个效应量在当时的 n 下只有 30% 的把握能检出来。从那以后我给自己定了条规矩任何 GLRT 结论必须附带一组功效曲线或者至少说明最小可检测效应量。这个习惯看起来多花十分钟但能避免你对着一个没有分辨能力的统计量瞎解释。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取