简介面向统计信号处理学习者与科研人员的广义最大似然比检验GLRTMATLAB仿真资源聚焦弱信号检测与噪声背景下异常判断问题适合正在学习假设检验、需要动手验证理论的本科高年级或研究生。压缩包仅3KB包含7个功能明确的m脚本分别覆盖仿真数据生成、似然函数与似然比计算、临界值求解、检测概率与虚警率评估等关键环节并基于连续波信号模型进行设置便于对照理论逐步复现。已有452人学习代码结构清晰、注释简明可帮助读者快速理解GLRT的决策流程观察不同参数对检测性能的影响并为实际雷达或通信场景中的检测问题提供可修改的仿真框架。同时资源围绕一个典型算例完整呈现从建立假设、计算似然比到统计判决的整个过程适合作为课程设计或论文复现的参考。1. 看到这串文件名先别急着解压MLE、GLRT、似然比检验说的是同一件事看到这串文件名做统计信号处理的人多半会心一笑。MLE、GLRT、Statistical test、似然比检验——四个词指向同一个技术方向用最大似然估计构造广义似然比检验回答信号到底存不存在、模型该不该换。这个方向的实际价值很直接只要观测能被某个概率分布描述想在两个假设之间做选择似然比检验就是最通用的框架之一参数未知时GLRT先用MLE把未知量估出来再比代价是渐近分布只在样本量足够时才准确。这篇笔记按原理→代码→参数→踩坑→验证展开适合正在做信号检测、模型对比或需要证明差异显著的从业者。新手跟第3章能跑通最小例子熟手直接看第4、5章的参数与坑位。2. 似然比检验的原理MLE是地基GLRT是处理未知参数的升级版2.1 最大似然估计为什么似然比检验离不开它回到最原始的问题。拿到一组观测 $x_1, x_2, \dots, x_n$假设它们独立同分布服从带参数 $\theta$ 的分布 $f(x;\theta)$。似然函数写成$$L(\theta) \prod_{i1}^{n} f(x_i; \theta)$$最大似然估计就是找出让 $L(\theta)$ 最大的那个 $\hat{\theta}$。这个概念看起来只是在做参数估计但它同时给统计假设检验提供了一把尺子如果某个假设能让数据出现的概率最大那这个假设就是相对合理的。于是检验问题被转化成比较受约束的假设和不受约束的假设各自对数据的解释能力而解释能力的度量就是似然函数值本身。似然比检验的统计量写成$$\Lambda(x) \frac{\max_{\theta \in \Theta_0} L(\theta)}{\max_{\theta \in \Theta} L(\theta)}$$分子是零假设约束下的最大似然值分母是全局最大似然值。因为 $\Theta_0$ 是 $\Theta$ 的子集分子永远不超过分母所以 $\Lambda$ 落在 $[0,1]$ 区间内。$\Lambda$ 越接近1说明零假设下数据也被解释得不错倾向于接受 $H_0$$\Lambda$ 明显小于1说明零假设太勉强应该拒绝。这里嵌套模型是整个方法的基石——如果两个假设不是嵌套关系$\Lambda$ 可能大于1整个框架立刻失效。实际计算时没人直接用 $\Lambda$原因在 Wilks 定理里。这个定理说的是在 $H_0$ 成立、样本量趋于无穷时$-2\ln\Lambda$ 渐近服从自由度为 $d$ 的卡方分布其中 $d \dim(\Theta_1) - \dim(\Theta_0)$。它的工程价值极大你不需要为每个具体问题做蒙特卡洛仿真找门限直接用卡方分布的分位数就能判决用scipy.stats.chi2.ppf和sf就能拿到门限和 p 值。代价是它是渐近结论样本量不够大、参数落在边界、模型不满足正则条件时卡方近似都会偏离这留给第5章细说。2.2 GLRT的定义把未知参数替换成MLE之后统计量发生了什么变化Neyman-Pearson 引理给出了最优检验前提是两个假设下的概率密度完全已知包括所有参数值。现实里这个前提几乎不成立——你多半只大概知道噪声方差信号幅度未知或者方差也要估计。广义似然比检验的做法很直白未知参数就用 MLE 估出来然后当作已知参数代入似然比$$t(x) \frac{f(x; \hat{\theta}_1)}{f(x; \hat{\theta}_0)}$$注意 $\hat{\theta}_0$ 是受 $H_0$ 约束的 MLE$\hat{\theta}_1$ 是无约束的 MLE两者不是一回事别在代码里省掉其中一个。直觉上因为估计过程偷看了数据GLRT 会比理想似然比检验更激进倾向于拒绝 $H_0$。这个偏差在渐近意义下会消失这是 Wilks 定理保证的但有限样本下必须靠仿真校准。以检测高斯白噪声中的直流信号为例接收模型是 $x_i \mu w_i$其中 $w_i \sim \mathcal{N}(0, \sigma^2)$。检验目标是 $H_0: \mu 0$ 对 $H_1: \mu eq 0$。当方差已知时$\mu$ 的 MLE 就是样本均值 $\bar{x}$这是高斯似然求导后得到的闭式解不需要迭代优化。把 MLE 代入似然比并整理$-2\ln\Lambda$ 化简成$$t(x) \frac{n\bar{x}^2}{\sigma^2}$$这个形式有清晰的物理含义样本均值偏离零假设的程度除以噪声方差再乘上样本量。偏离越大、数据量越多越怀疑信号存在。它还说明了一件事——在高斯加性噪声模型下GLRT 最终等价于能量检测器或均值检测器统计量是样本均值的函数不是更复杂的东西。如果只关心信号方向比如只检测正的直流偏移就把判决改成单边不仅要 $t(x)$ 超过门限还要求 $\bar{x} 0$。2.3 自由度与模型嵌套这个细节决定了p值对不对自由度 $d$ 是两个参数空间维度之差不是被检验参数的个数。检验两个独立样本组的均值是否相等时$H_1$ 的参数空间是 $(\mu_1, \mu_2, \sigma^2)$维度3$H_0$ 下 $\mu_1 \mu_2 \mu$参数空间是 $(\mu, \sigma^2)$维度2自由度是 $3-21$。如果你凭直觉认为两个均值就是两个参数自由度应该是2那 p 值就会系统性偏大检验变得过于保守本来显著的差异可能被判成不显著。嵌套性是另一个高频出错点。嵌套的意思是 $H_0$ 的参数空间是 $H_1$ 参数空间的子集。检验 $\mu0$ 对 $\mu eq 0$ 是嵌套的检验观测服从正态分布对观测服从指数分布就不是嵌套的两者不能做似然比检验。非嵌套模型在工程上更常用的比较工具是 AIC、BIC 这类信息准则它们的推导出发点不同不能混用。还有一类边界问题比如 $H_0: \mu \ge 0$ 而真实参数恰好落在边界 $\mu0$ 上此时卡方近似的收敛速度会明显变慢实际虚警率比名义水平偏大这在第5章的 5.4 节有对应的处理方法。3. 用GLRT检测高斯噪声中的直流信号最小可运行的Python实现3.1 先写清楚假设再写代码写统计检验代码前第一件事是用注释把三件事写清楚数据模型、零假设与备择假设、哪些参数已知。这个习惯能省掉大量返工。本文的场景如下数据模型是 $x_i \mu w_i$$i 1, \dots, n$$w_i$ 独立同分布服从 $\mathcal{N}(0, \sigma^2)$检验目标是 $H_0: \mu 0$ 对 $H_1: \mu eq 0$先假定方差 $\sigma^2$ 已知。为什么先假定方差已知因为这能让 GLRT 统计量有解析形式并且在 $H_0$ 下精确服从 $\chi^2_1$ 分布。你可以用蒙特卡洛仿真验证代码写对了没有。如果一上来就假定方差未知统计量的精确分布变成 $F(1, n-1)$代码跑出的结果对不对就很难判断——到底是理论错了还是实现错了先建立一个可信基线再逐步放开假设这是工程上更稳的推进方式。实际计算时$H_1$ 下 $\mu$ 的 MLE 就是样本均值 $\bar{x} \frac{1}{n}\sum_{i1}^{n}x_i$这是高斯似然对 $\mu$ 求导得到的闭式解。有人会用scipy.optimize.minimize数值最大化似然函数结果和np.mean(x)一致但慢了几十倍。高斯模型这类可解析求解的场景直接用闭式解节省时间也避免数值问题。3.2 核心代码GLRT统计量、p值与判决import numpy as np from scipy import stats def glrt_dc_detect(x, sigma2, alpha0.05): 高斯白噪声中检测未知直流信号的GLRT。 参数 ---- x : array_like 观测序列 sigma2 : float 噪声方差(已知) alpha : float 显著性水平即允许的虚警概率 返回 ---- t_stat : float GLRT统计量 n * mean^2 / sigma2 p_value : float 零假设下出现当前或更极端统计量的概率 reject : bool True 表示拒绝 H0认为存在直流信号 x np.asarray(x, dtypefloat) n x.size # H1下 mu 的MLE高斯似然的闭式解是样本均值 mu_hat np.mean(x) # GLRT统计量-2*ln(L0/L1) 化简后的解析形式 # 不直接计算两个似然函数再相除避开数值下溢 t_stat n * mu_hat**2 / sigma2 # 卡方分布自由度1双边检验 mu ! 0 p_value stats.chi2.sf(t_stat, df1) threshold stats.chi2.ppf(1 - alpha, df1) return t_stat, p_value, t_stat threshold这里的三个设计选择值得说明。第一统计量用化简式而不是原始似然比定义当 $n$ 超过几十直接用prod(f(x_i))计算似然函数乘积会下溢成0取对数再相减同样会损失精度化简式完全绕开这个问题。第二chi2.sf算右尾概率正好对应 p 值的定义——零假设下出现当前或更极端统计量的概率门限用ppf(1-alpha, df1)取上侧分位数两侧逻辑一致。第三函数返回统计量、p 值、判决结果三个值后续做 Monte Carlo 仿真时直接取用不需要重复计算。3.3 Monte Carlo验证5千次试验看虚警率验证代码正确性的最好方法不是反复核对理论推导而是做仿真。思路很简单在 $H_0$ 下生成大量数据集每次调用glrt_dc_detect统计误报比例看它是否接近 $\alpha$。def estimate_false_alarm(n50, sigma21.0, alpha0.05, n_trials5000, seed42): rng np.random.default_rng(seed) false_alarms 0 for _ in range(n_trials): x rng.normal(0.0, np.sqrt(sigma2), n) _, _, reject glrt_dc_detect(x, sigma2, alpha) false_alarms int(reject) return false_alarms / n_trials for a in [0.01, 0.05, 0.1]: fa estimate_false_alarm(alphaa) print(falpha{a:.2f}, 经验虚警率{fa:.4f})几个参数值得细说。n_trials5000是最低可接受的仿真次数经验频率的标准误差约为 $\sqrt{p(1-p)/N}$在 $p0.05$ 时约0.003足够判断虚警率是否显著偏离理论值。seed42不是随便选的固定随机种子才能让结果可复现同事跑出来的数字和你的一致。如果两次运行结果不同不是代码有 bug而是没有固定种子。rng.normal走 NumPy 1.17 之后推荐的default_rng接口不要再用旧的全局np.random.seed加np.random.normal混着写。运行后你会看到经验虚警率在理论值附近小幅波动。如果偏差超过0.01第一步检查自由度是否写错第二步检查统计量是不是忘了乘 $n$第三步检查H0下的数据生成有没有混入非零均值。这三个低级错误几乎覆盖了所有虚警率对不上的场景。提示把这段虚警率自检代码保留在你的工具库里。以后每次改动统计量、换数据生成方式、调自由度都先跑一遍。这是检验代码回归测试的核心比任何代码审查都可靠。3.4 检测概率仿真扫描信噪比画性能曲线虚警率验证通过只说明零假设下没有失控还要看备择假设下能不能检测出来。固定虚警率检测概率与信噪比之间的关系就是这个检验的性能标尺。def monte_carlo_pd(snr_db, n50, sigma21.0, alpha0.05, n_trials5000, seed0): 给定信噪比(dB)下的检测概率。 snr_db 10*log10(mu^2 / sigma2) rng np.random.default_rng(seed) mu np.sqrt(sigma2 * 10**(snr_db / 10.0)) detects 0 for _ in range(n_trials): x rng.normal(mu, np.sqrt(sigma2), n) _, _, reject glrt_dc_detect(x, sigma2, alphaalpha) detects int(reject) return detects / n_trials for snr in [-10, -5, 0, 5, 10]: pd monte_carlo_pd(snr) print(fSNR{snr:3d} dB, P_d{pd:.3f})信噪比定义成 $10\log_{10}(\mu^2/\sigma^2)$单位是 dB。负10dB 时信号功率是噪声功率的十分之一检测概率应该接近 $\alpha$基本靠猜0dB 时均值幅度等于噪声标准差检测概率明显上升10dB 以上几乎必然检测到。如果曲线不符合这个趋势多半是数据生成或统计量构造有错。常见错误是把 $\mu$ 生成成 $\sqrt{\sigma^2 \cdot 10^{snr/10}}$ 时忘记外层开根号或者把 $\sigma^2$ 错当 $\sigma$ 用。这个函数本身也可复用换场景时只需要改数据生成那一行性能评估框架不用动。4. 参数怎么设显著性水平、自由度与样本量的工程选择4.1 显著性水平不要照抄0.05从应用场景倒推门限教科书里 0.05 几乎成了默认值实际工程里这个数值要由应用承担的成本决定。学术论文里 0.05 用于控制错误发现率尚可接受雷达检测、故障告警、医学诊断这类场景一次虚警的代价很高虚警率指标可能要求 $10^{-5}$ 甚至更低。此时用卡方分布理论分位数作为门限是否还可靠只要统计量在 $H_0$ 下确实服从卡方分布分位数本身就是精确的不会因为显著性水平变小而失效。真正的风险来自样本量不足尾部偏差在小 $\alpha$ 下会被放大小样本加小 $\alpha$ 的场景必须用经验门限。经验门限的做法是在 $H_0$ 下生成 $N$ 个数据集算出 $N$ 个统计量 $t_1, \dots, t_N$排序后取 $(1-\alpha)$ 分位数作为门限。这个门限不依赖任何渐近近似代价是需要几千到几万次仿真。我通常把理论卡方门限和经验门限一起打出来对比若偏差超过20%就采用经验门限。对需要反复使用的检验场景把经验门限预计算好存下来运行时查表性能完全可接受。4.2 自由度怎么数才不出错维度差方法自由度计算的唯一可靠方法是维度差$\dim(\Theta_1) - \dim(\Theta_0)$。以下是四种常见模型的具体数值检验场景H1参数空间与维度H0参数空间与维度自由度直流信号 mu0方差已知mu维度1无参数维度01直流信号 mu0方差未知mu, sigma2维度2sigma2维度11两组均值相等方差未知mu1, mu2, sigma2维度3mu, sigma2维度21线性回归全模型 vs 去掉2个系数原模型 p1 个参数减元模型 p-1 个参数2注意第二行和第三行的自由度都是1。方差未知时不要以为参数变多了所以自由度变大——两个假设下方差都要估计这个维度被抵消了。自由度真正变化的场景是 $H_1$ 比 $H_0$ 多估计了若干个系数比如回归模型里比较包含两个额外自变量与不包含的情况自由度才是2。写代码前把两边的参数空间列出来把维度差的数字写在函数注释里这是成本最低的防错手段。4.3 样本量与精确分布何时不能靠卡方近似Wilks 定理是渐近结果样本量多少才算足够大没有严格边界。以方差未知的直流检测为例GLRT 统计量在 $H_0$ 下精确服从 $F(1, n-1)$ 分布。$n20$ 时 $F(1,19)$ 的0.95分位数约为4.38而 $\chi^2_1$ 的0.95分位数是3.84——用卡方门限意味着实际显著性水平高于名义值虚警偏多。$n50$ 时约4.03差距缩到5%$n100$ 时约3.94差距约2.6%。工程上我的习惯是$n \ge 100$ 才放心用卡方近似$n$ 在20到100之间尽量用精确分布查scipy.stats.f.ppf(1-alpha, 1, n-1)做门限$n 20$ 时 GLRT 的小样本性质不稳定考虑置换检验这类重抽样方法。置换检验对分布假设要求低但每次检验需要几千次重抽样计算量大适合离线分析不适合在线实时检测。这里不展开实现细节但要记住它是一个可靠的备选项。4.4 模型假设的适用前提写代码前就要想清楚最后一个参数类问题不是数值而是模型设定。GLRT 适用前提包括样本独立同分布或至少似然函数可写参数空间光滑真实参数不在边界上模型嵌套。如果你的数据是时间序列且强相关或者备择假设的参数空间不光滑再调参数也救不回来。正确的做法是换检验框架比如基于秩的非参数检验或专门处理序列相关性的似然比修正。判断一个检验方法是否适用优先级高于调整参数——参数调得再好模型设错了也是白做。5. 避坑似然比检验最常见的5个翻车现场5.1 似然比大于1取对数直接报错现象代码算出 $\Lambda 1$np.log得到正数p 值变成负的或毫无意义。原因几乎只有两种一是两个模型不是嵌套关系分母不是全局最大化二是代码里分子分母写反了。解决先确认 $H_0$ 的参数空间确实是 $H_1$ 的子集再检查L0和L1的赋值顺序。我习惯在函数 docstring 里把分子分母的含义写死并在关键位置加一行assert lambda_ 1 1e-9一旦违反立刻暴露问题不等到后面算出奇怪数字才回头查。5.2 p值整体偏移经验虚警率对不上显著性水平现象$H_0$ 下仿真经验虚警率稳定在0.08而不是0.05。原因自由度写错、统计量公式漏了 $n$、或用了错误的方差值。排查顺序按三步走。第一步打印统计量的均值$H_0$ 下 $\chi^2_1$ 的均值是1如果均值明显偏离统计量本身就有问题。第二步核对仿真代码里的数据生成是否真的满足 $H_0$。第三步检查chi2.sf(t, df)里df是不是维度差而非参数个数。按这个顺序排查基本十分钟内可以定位。我一直保留一个 $H_0$ 仿真的最小脚本任何检验代码改动后都先跑一遍让回归测试替我把关。5.3 对数似然全是-inf统计量变成NaN现象计算对数似然时出现RuntimeWarning: divide by zero输出是nan或无穷大。原因直接连乘密度函数再取对数样本量一大就下溢或者密度函数在某个点取到0对数变成负无穷。解决所有似然计算一律使用对数似然并逐样本累加禁止np.prod。对高斯模型直接用scipy.stats.norm.logpdf(x, loc, scale).sum()由库函数处理尾部细节。如果必须自己写密度函数用np.log时要对参数范围做保护比如方差下限加一个eps1e-12。这个坑在做非高斯模型时尤其常见指数族之外的概率密度很容易在某个点取0。5.4 真实虚警率比理论值偏大重复多次都一样现象$H_0$ 下经验虚警率稳定高于 $\alpha$并且和样本量关系不大。原因参数位于参数空间边界。典型场景是单边检验 $H_0: \mu \ge 0$ 里真实 $\mu0$ 正好落在边界上或方差检验中 $H_0: \sigma^20$。边界破坏了 Wilks 定理所需的正则条件卡方近似不成立。解决改用仿真校准门限或者用专门处理边界问题的混合卡方分布。工程上如果只是要一个判决结论做 $10^5$ 次仿真标定经验门限代码简单结论可靠。别试图用更大的样本量硬撑——在某些边界场景下样本量再大也救不回卡方近似。5.5 换了随机种子检测概率结果差很多现象seed0时检测概率 $P_d0.65$seed1时 $P_d0.71$不知道该信谁。原因仿真次数太少估计标准差太大。$P_d \approx 0.68$ 时1000 次试验的标准差约0.015但统计量接近门限时单次试验结果对种子更敏感实际抖动会更大。解决试验次数提到5000以上固定种子并把不同种子下的结果都打出来做一致性检查应该在 ±0.02 范围内波动。如果项目对性能曲线精度要求高额外保存原始统计量列表而不是只存一个均值。固定种子不是为了好看是为了让结论可复现——否则同事复现你的结果时会对不上最后只能花时间排查一个不存在的问题。6. 收尾与你的实际场景对接以及三个必须做的自检6.1 拿到类似的项目包先看什么如果你手上的压缩包是统计检验代码的常见形态先别急着改代码按三个层次读。第一层跑通 demo 脚本复现它输出的数字第二层看核心函数签名返回什么、假定哪些参数已知判断能不能直接用你的数据喂进去第三层找到数据生成与假设设定的注释确认与你的场景一致。spellcw5这类后缀通常是作者或版本的标识不影响使用。按正常流程解压即可代码包的完整性和安全性检查是另一件事这里不展开。6.2 三个自检方法改动任何假设后都跑一遍第一个自检$H_0$ 下经验虚警率与名义 $\alpha$ 一致。第二个自检$H_1$ 下固定信噪比检测概率随样本量 $n$ 增大单调上升。第三个自检统计量的经验直方图与理论卡方密度叠加后形状吻合直观确认渐近逼近成立。这三个自检加起来不到50行能拦下绝大多数翻车。我自己的习惯是把自检写成独立脚本不混进业务代码每次改完假设和参数单独重跑。最后分享一个长期习惯每写完一组检验代码把数据生成、假设设定、自由度计算三件事单独抽成模板。下次换场景先回答三个问题——方差已知还是未知、单边还是双边、模型是否嵌套——回答完这三个问题检验代码的框架基本不用动。这套流程帮我处理过几十个检测与对比任务多数统计错误都源于这三个问题没想清就动手写。希望帮到你。本文还有配套的精品资源点击获取 最大似然估计与广义似然比检验：从原理到Python工程实践
简介面向统计信号处理学习者与科研人员的广义最大似然比检验GLRTMATLAB仿真资源聚焦弱信号检测与噪声背景下异常判断问题适合正在学习假设检验、需要动手验证理论的本科高年级或研究生。压缩包仅3KB包含7个功能明确的m脚本分别覆盖仿真数据生成、似然函数与似然比计算、临界值求解、检测概率与虚警率评估等关键环节并基于连续波信号模型进行设置便于对照理论逐步复现。已有452人学习代码结构清晰、注释简明可帮助读者快速理解GLRT的决策流程观察不同参数对检测性能的影响并为实际雷达或通信场景中的检测问题提供可修改的仿真框架。同时资源围绕一个典型算例完整呈现从建立假设、计算似然比到统计判决的整个过程适合作为课程设计或论文复现的参考。1. 看到这串文件名先别急着解压MLE、GLRT、似然比检验说的是同一件事看到这串文件名做统计信号处理的人多半会心一笑。MLE、GLRT、Statistical test、似然比检验——四个词指向同一个技术方向用最大似然估计构造广义似然比检验回答信号到底存不存在、模型该不该换。这个方向的实际价值很直接只要观测能被某个概率分布描述想在两个假设之间做选择似然比检验就是最通用的框架之一参数未知时GLRT先用MLE把未知量估出来再比代价是渐近分布只在样本量足够时才准确。这篇笔记按原理→代码→参数→踩坑→验证展开适合正在做信号检测、模型对比或需要证明差异显著的从业者。新手跟第3章能跑通最小例子熟手直接看第4、5章的参数与坑位。2. 似然比检验的原理MLE是地基GLRT是处理未知参数的升级版2.1 最大似然估计为什么似然比检验离不开它回到最原始的问题。拿到一组观测 $x_1, x_2, \dots, x_n$假设它们独立同分布服从带参数 $\theta$ 的分布 $f(x;\theta)$。似然函数写成$$L(\theta) \prod_{i1}^{n} f(x_i; \theta)$$最大似然估计就是找出让 $L(\theta)$ 最大的那个 $\hat{\theta}$。这个概念看起来只是在做参数估计但它同时给统计假设检验提供了一把尺子如果某个假设能让数据出现的概率最大那这个假设就是相对合理的。于是检验问题被转化成比较受约束的假设和不受约束的假设各自对数据的解释能力而解释能力的度量就是似然函数值本身。似然比检验的统计量写成$$\Lambda(x) \frac{\max_{\theta \in \Theta_0} L(\theta)}{\max_{\theta \in \Theta} L(\theta)}$$分子是零假设约束下的最大似然值分母是全局最大似然值。因为 $\Theta_0$ 是 $\Theta$ 的子集分子永远不超过分母所以 $\Lambda$ 落在 $[0,1]$ 区间内。$\Lambda$ 越接近1说明零假设下数据也被解释得不错倾向于接受 $H_0$$\Lambda$ 明显小于1说明零假设太勉强应该拒绝。这里嵌套模型是整个方法的基石——如果两个假设不是嵌套关系$\Lambda$ 可能大于1整个框架立刻失效。实际计算时没人直接用 $\Lambda$原因在 Wilks 定理里。这个定理说的是在 $H_0$ 成立、样本量趋于无穷时$-2\ln\Lambda$ 渐近服从自由度为 $d$ 的卡方分布其中 $d \dim(\Theta_1) - \dim(\Theta_0)$。它的工程价值极大你不需要为每个具体问题做蒙特卡洛仿真找门限直接用卡方分布的分位数就能判决用scipy.stats.chi2.ppf和sf就能拿到门限和 p 值。代价是它是渐近结论样本量不够大、参数落在边界、模型不满足正则条件时卡方近似都会偏离这留给第5章细说。2.2 GLRT的定义把未知参数替换成MLE之后统计量发生了什么变化Neyman-Pearson 引理给出了最优检验前提是两个假设下的概率密度完全已知包括所有参数值。现实里这个前提几乎不成立——你多半只大概知道噪声方差信号幅度未知或者方差也要估计。广义似然比检验的做法很直白未知参数就用 MLE 估出来然后当作已知参数代入似然比$$t(x) \frac{f(x; \hat{\theta}_1)}{f(x; \hat{\theta}_0)}$$注意 $\hat{\theta}_0$ 是受 $H_0$ 约束的 MLE$\hat{\theta}_1$ 是无约束的 MLE两者不是一回事别在代码里省掉其中一个。直觉上因为估计过程偷看了数据GLRT 会比理想似然比检验更激进倾向于拒绝 $H_0$。这个偏差在渐近意义下会消失这是 Wilks 定理保证的但有限样本下必须靠仿真校准。以检测高斯白噪声中的直流信号为例接收模型是 $x_i \mu w_i$其中 $w_i \sim \mathcal{N}(0, \sigma^2)$。检验目标是 $H_0: \mu 0$ 对 $H_1: \mu eq 0$。当方差已知时$\mu$ 的 MLE 就是样本均值 $\bar{x}$这是高斯似然求导后得到的闭式解不需要迭代优化。把 MLE 代入似然比并整理$-2\ln\Lambda$ 化简成$$t(x) \frac{n\bar{x}^2}{\sigma^2}$$这个形式有清晰的物理含义样本均值偏离零假设的程度除以噪声方差再乘上样本量。偏离越大、数据量越多越怀疑信号存在。它还说明了一件事——在高斯加性噪声模型下GLRT 最终等价于能量检测器或均值检测器统计量是样本均值的函数不是更复杂的东西。如果只关心信号方向比如只检测正的直流偏移就把判决改成单边不仅要 $t(x)$ 超过门限还要求 $\bar{x} 0$。2.3 自由度与模型嵌套这个细节决定了p值对不对自由度 $d$ 是两个参数空间维度之差不是被检验参数的个数。检验两个独立样本组的均值是否相等时$H_1$ 的参数空间是 $(\mu_1, \mu_2, \sigma^2)$维度3$H_0$ 下 $\mu_1 \mu_2 \mu$参数空间是 $(\mu, \sigma^2)$维度2自由度是 $3-21$。如果你凭直觉认为两个均值就是两个参数自由度应该是2那 p 值就会系统性偏大检验变得过于保守本来显著的差异可能被判成不显著。嵌套性是另一个高频出错点。嵌套的意思是 $H_0$ 的参数空间是 $H_1$ 参数空间的子集。检验 $\mu0$ 对 $\mu eq 0$ 是嵌套的检验观测服从正态分布对观测服从指数分布就不是嵌套的两者不能做似然比检验。非嵌套模型在工程上更常用的比较工具是 AIC、BIC 这类信息准则它们的推导出发点不同不能混用。还有一类边界问题比如 $H_0: \mu \ge 0$ 而真实参数恰好落在边界 $\mu0$ 上此时卡方近似的收敛速度会明显变慢实际虚警率比名义水平偏大这在第5章的 5.4 节有对应的处理方法。3. 用GLRT检测高斯噪声中的直流信号最小可运行的Python实现3.1 先写清楚假设再写代码写统计检验代码前第一件事是用注释把三件事写清楚数据模型、零假设与备择假设、哪些参数已知。这个习惯能省掉大量返工。本文的场景如下数据模型是 $x_i \mu w_i$$i 1, \dots, n$$w_i$ 独立同分布服从 $\mathcal{N}(0, \sigma^2)$检验目标是 $H_0: \mu 0$ 对 $H_1: \mu eq 0$先假定方差 $\sigma^2$ 已知。为什么先假定方差已知因为这能让 GLRT 统计量有解析形式并且在 $H_0$ 下精确服从 $\chi^2_1$ 分布。你可以用蒙特卡洛仿真验证代码写对了没有。如果一上来就假定方差未知统计量的精确分布变成 $F(1, n-1)$代码跑出的结果对不对就很难判断——到底是理论错了还是实现错了先建立一个可信基线再逐步放开假设这是工程上更稳的推进方式。实际计算时$H_1$ 下 $\mu$ 的 MLE 就是样本均值 $\bar{x} \frac{1}{n}\sum_{i1}^{n}x_i$这是高斯似然对 $\mu$ 求导得到的闭式解。有人会用scipy.optimize.minimize数值最大化似然函数结果和np.mean(x)一致但慢了几十倍。高斯模型这类可解析求解的场景直接用闭式解节省时间也避免数值问题。3.2 核心代码GLRT统计量、p值与判决import numpy as np from scipy import stats def glrt_dc_detect(x, sigma2, alpha0.05): 高斯白噪声中检测未知直流信号的GLRT。 参数 ---- x : array_like 观测序列 sigma2 : float 噪声方差(已知) alpha : float 显著性水平即允许的虚警概率 返回 ---- t_stat : float GLRT统计量 n * mean^2 / sigma2 p_value : float 零假设下出现当前或更极端统计量的概率 reject : bool True 表示拒绝 H0认为存在直流信号 x np.asarray(x, dtypefloat) n x.size # H1下 mu 的MLE高斯似然的闭式解是样本均值 mu_hat np.mean(x) # GLRT统计量-2*ln(L0/L1) 化简后的解析形式 # 不直接计算两个似然函数再相除避开数值下溢 t_stat n * mu_hat**2 / sigma2 # 卡方分布自由度1双边检验 mu ! 0 p_value stats.chi2.sf(t_stat, df1) threshold stats.chi2.ppf(1 - alpha, df1) return t_stat, p_value, t_stat threshold这里的三个设计选择值得说明。第一统计量用化简式而不是原始似然比定义当 $n$ 超过几十直接用prod(f(x_i))计算似然函数乘积会下溢成0取对数再相减同样会损失精度化简式完全绕开这个问题。第二chi2.sf算右尾概率正好对应 p 值的定义——零假设下出现当前或更极端统计量的概率门限用ppf(1-alpha, df1)取上侧分位数两侧逻辑一致。第三函数返回统计量、p 值、判决结果三个值后续做 Monte Carlo 仿真时直接取用不需要重复计算。3.3 Monte Carlo验证5千次试验看虚警率验证代码正确性的最好方法不是反复核对理论推导而是做仿真。思路很简单在 $H_0$ 下生成大量数据集每次调用glrt_dc_detect统计误报比例看它是否接近 $\alpha$。def estimate_false_alarm(n50, sigma21.0, alpha0.05, n_trials5000, seed42): rng np.random.default_rng(seed) false_alarms 0 for _ in range(n_trials): x rng.normal(0.0, np.sqrt(sigma2), n) _, _, reject glrt_dc_detect(x, sigma2, alpha) false_alarms int(reject) return false_alarms / n_trials for a in [0.01, 0.05, 0.1]: fa estimate_false_alarm(alphaa) print(falpha{a:.2f}, 经验虚警率{fa:.4f})几个参数值得细说。n_trials5000是最低可接受的仿真次数经验频率的标准误差约为 $\sqrt{p(1-p)/N}$在 $p0.05$ 时约0.003足够判断虚警率是否显著偏离理论值。seed42不是随便选的固定随机种子才能让结果可复现同事跑出来的数字和你的一致。如果两次运行结果不同不是代码有 bug而是没有固定种子。rng.normal走 NumPy 1.17 之后推荐的default_rng接口不要再用旧的全局np.random.seed加np.random.normal混着写。运行后你会看到经验虚警率在理论值附近小幅波动。如果偏差超过0.01第一步检查自由度是否写错第二步检查统计量是不是忘了乘 $n$第三步检查H0下的数据生成有没有混入非零均值。这三个低级错误几乎覆盖了所有虚警率对不上的场景。提示把这段虚警率自检代码保留在你的工具库里。以后每次改动统计量、换数据生成方式、调自由度都先跑一遍。这是检验代码回归测试的核心比任何代码审查都可靠。3.4 检测概率仿真扫描信噪比画性能曲线虚警率验证通过只说明零假设下没有失控还要看备择假设下能不能检测出来。固定虚警率检测概率与信噪比之间的关系就是这个检验的性能标尺。def monte_carlo_pd(snr_db, n50, sigma21.0, alpha0.05, n_trials5000, seed0): 给定信噪比(dB)下的检测概率。 snr_db 10*log10(mu^2 / sigma2) rng np.random.default_rng(seed) mu np.sqrt(sigma2 * 10**(snr_db / 10.0)) detects 0 for _ in range(n_trials): x rng.normal(mu, np.sqrt(sigma2), n) _, _, reject glrt_dc_detect(x, sigma2, alphaalpha) detects int(reject) return detects / n_trials for snr in [-10, -5, 0, 5, 10]: pd monte_carlo_pd(snr) print(fSNR{snr:3d} dB, P_d{pd:.3f})信噪比定义成 $10\log_{10}(\mu^2/\sigma^2)$单位是 dB。负10dB 时信号功率是噪声功率的十分之一检测概率应该接近 $\alpha$基本靠猜0dB 时均值幅度等于噪声标准差检测概率明显上升10dB 以上几乎必然检测到。如果曲线不符合这个趋势多半是数据生成或统计量构造有错。常见错误是把 $\mu$ 生成成 $\sqrt{\sigma^2 \cdot 10^{snr/10}}$ 时忘记外层开根号或者把 $\sigma^2$ 错当 $\sigma$ 用。这个函数本身也可复用换场景时只需要改数据生成那一行性能评估框架不用动。4. 参数怎么设显著性水平、自由度与样本量的工程选择4.1 显著性水平不要照抄0.05从应用场景倒推门限教科书里 0.05 几乎成了默认值实际工程里这个数值要由应用承担的成本决定。学术论文里 0.05 用于控制错误发现率尚可接受雷达检测、故障告警、医学诊断这类场景一次虚警的代价很高虚警率指标可能要求 $10^{-5}$ 甚至更低。此时用卡方分布理论分位数作为门限是否还可靠只要统计量在 $H_0$ 下确实服从卡方分布分位数本身就是精确的不会因为显著性水平变小而失效。真正的风险来自样本量不足尾部偏差在小 $\alpha$ 下会被放大小样本加小 $\alpha$ 的场景必须用经验门限。经验门限的做法是在 $H_0$ 下生成 $N$ 个数据集算出 $N$ 个统计量 $t_1, \dots, t_N$排序后取 $(1-\alpha)$ 分位数作为门限。这个门限不依赖任何渐近近似代价是需要几千到几万次仿真。我通常把理论卡方门限和经验门限一起打出来对比若偏差超过20%就采用经验门限。对需要反复使用的检验场景把经验门限预计算好存下来运行时查表性能完全可接受。4.2 自由度怎么数才不出错维度差方法自由度计算的唯一可靠方法是维度差$\dim(\Theta_1) - \dim(\Theta_0)$。以下是四种常见模型的具体数值检验场景H1参数空间与维度H0参数空间与维度自由度直流信号 mu0方差已知mu维度1无参数维度01直流信号 mu0方差未知mu, sigma2维度2sigma2维度11两组均值相等方差未知mu1, mu2, sigma2维度3mu, sigma2维度21线性回归全模型 vs 去掉2个系数原模型 p1 个参数减元模型 p-1 个参数2注意第二行和第三行的自由度都是1。方差未知时不要以为参数变多了所以自由度变大——两个假设下方差都要估计这个维度被抵消了。自由度真正变化的场景是 $H_1$ 比 $H_0$ 多估计了若干个系数比如回归模型里比较包含两个额外自变量与不包含的情况自由度才是2。写代码前把两边的参数空间列出来把维度差的数字写在函数注释里这是成本最低的防错手段。4.3 样本量与精确分布何时不能靠卡方近似Wilks 定理是渐近结果样本量多少才算足够大没有严格边界。以方差未知的直流检测为例GLRT 统计量在 $H_0$ 下精确服从 $F(1, n-1)$ 分布。$n20$ 时 $F(1,19)$ 的0.95分位数约为4.38而 $\chi^2_1$ 的0.95分位数是3.84——用卡方门限意味着实际显著性水平高于名义值虚警偏多。$n50$ 时约4.03差距缩到5%$n100$ 时约3.94差距约2.6%。工程上我的习惯是$n \ge 100$ 才放心用卡方近似$n$ 在20到100之间尽量用精确分布查scipy.stats.f.ppf(1-alpha, 1, n-1)做门限$n 20$ 时 GLRT 的小样本性质不稳定考虑置换检验这类重抽样方法。置换检验对分布假设要求低但每次检验需要几千次重抽样计算量大适合离线分析不适合在线实时检测。这里不展开实现细节但要记住它是一个可靠的备选项。4.4 模型假设的适用前提写代码前就要想清楚最后一个参数类问题不是数值而是模型设定。GLRT 适用前提包括样本独立同分布或至少似然函数可写参数空间光滑真实参数不在边界上模型嵌套。如果你的数据是时间序列且强相关或者备择假设的参数空间不光滑再调参数也救不回来。正确的做法是换检验框架比如基于秩的非参数检验或专门处理序列相关性的似然比修正。判断一个检验方法是否适用优先级高于调整参数——参数调得再好模型设错了也是白做。5. 避坑似然比检验最常见的5个翻车现场5.1 似然比大于1取对数直接报错现象代码算出 $\Lambda 1$np.log得到正数p 值变成负的或毫无意义。原因几乎只有两种一是两个模型不是嵌套关系分母不是全局最大化二是代码里分子分母写反了。解决先确认 $H_0$ 的参数空间确实是 $H_1$ 的子集再检查L0和L1的赋值顺序。我习惯在函数 docstring 里把分子分母的含义写死并在关键位置加一行assert lambda_ 1 1e-9一旦违反立刻暴露问题不等到后面算出奇怪数字才回头查。5.2 p值整体偏移经验虚警率对不上显著性水平现象$H_0$ 下仿真经验虚警率稳定在0.08而不是0.05。原因自由度写错、统计量公式漏了 $n$、或用了错误的方差值。排查顺序按三步走。第一步打印统计量的均值$H_0$ 下 $\chi^2_1$ 的均值是1如果均值明显偏离统计量本身就有问题。第二步核对仿真代码里的数据生成是否真的满足 $H_0$。第三步检查chi2.sf(t, df)里df是不是维度差而非参数个数。按这个顺序排查基本十分钟内可以定位。我一直保留一个 $H_0$ 仿真的最小脚本任何检验代码改动后都先跑一遍让回归测试替我把关。5.3 对数似然全是-inf统计量变成NaN现象计算对数似然时出现RuntimeWarning: divide by zero输出是nan或无穷大。原因直接连乘密度函数再取对数样本量一大就下溢或者密度函数在某个点取到0对数变成负无穷。解决所有似然计算一律使用对数似然并逐样本累加禁止np.prod。对高斯模型直接用scipy.stats.norm.logpdf(x, loc, scale).sum()由库函数处理尾部细节。如果必须自己写密度函数用np.log时要对参数范围做保护比如方差下限加一个eps1e-12。这个坑在做非高斯模型时尤其常见指数族之外的概率密度很容易在某个点取0。5.4 真实虚警率比理论值偏大重复多次都一样现象$H_0$ 下经验虚警率稳定高于 $\alpha$并且和样本量关系不大。原因参数位于参数空间边界。典型场景是单边检验 $H_0: \mu \ge 0$ 里真实 $\mu0$ 正好落在边界上或方差检验中 $H_0: \sigma^20$。边界破坏了 Wilks 定理所需的正则条件卡方近似不成立。解决改用仿真校准门限或者用专门处理边界问题的混合卡方分布。工程上如果只是要一个判决结论做 $10^5$ 次仿真标定经验门限代码简单结论可靠。别试图用更大的样本量硬撑——在某些边界场景下样本量再大也救不回卡方近似。5.5 换了随机种子检测概率结果差很多现象seed0时检测概率 $P_d0.65$seed1时 $P_d0.71$不知道该信谁。原因仿真次数太少估计标准差太大。$P_d \approx 0.68$ 时1000 次试验的标准差约0.015但统计量接近门限时单次试验结果对种子更敏感实际抖动会更大。解决试验次数提到5000以上固定种子并把不同种子下的结果都打出来做一致性检查应该在 ±0.02 范围内波动。如果项目对性能曲线精度要求高额外保存原始统计量列表而不是只存一个均值。固定种子不是为了好看是为了让结论可复现——否则同事复现你的结果时会对不上最后只能花时间排查一个不存在的问题。6. 收尾与你的实际场景对接以及三个必须做的自检6.1 拿到类似的项目包先看什么如果你手上的压缩包是统计检验代码的常见形态先别急着改代码按三个层次读。第一层跑通 demo 脚本复现它输出的数字第二层看核心函数签名返回什么、假定哪些参数已知判断能不能直接用你的数据喂进去第三层找到数据生成与假设设定的注释确认与你的场景一致。spellcw5这类后缀通常是作者或版本的标识不影响使用。按正常流程解压即可代码包的完整性和安全性检查是另一件事这里不展开。6.2 三个自检方法改动任何假设后都跑一遍第一个自检$H_0$ 下经验虚警率与名义 $\alpha$ 一致。第二个自检$H_1$ 下固定信噪比检测概率随样本量 $n$ 增大单调上升。第三个自检统计量的经验直方图与理论卡方密度叠加后形状吻合直观确认渐近逼近成立。这三个自检加起来不到50行能拦下绝大多数翻车。我自己的习惯是把自检写成独立脚本不混进业务代码每次改完假设和参数单独重跑。最后分享一个长期习惯每写完一组检验代码把数据生成、假设设定、自由度计算三件事单独抽成模板。下次换场景先回答三个问题——方差已知还是未知、单边还是双边、模型是否嵌套——回答完这三个问题检验代码的框架基本不用动。这套流程帮我处理过几十个检测与对比任务多数统计错误都源于这三个问题没想清就动手写。希望帮到你。本文还有配套的精品资源点击获取