简介这份资源是面向统计信号处理学习者与科研人员的MATLAB仿真代码包聚焦广义最大似然比检验GLRT这一假设检验方法帮助理解弱信号与噪声环境下的检测原理。压缩包共7个文件全部为m脚本整体约3KB涵盖数据生成、似然函数计算、临界值确定、决策规则与性能评估等模块可完整复现GLRT的仿真流程。已有452人学习下载适合希望从理论走向实践、掌握检测器设计的中高级读者。通过运行这些脚本读者能直观观察不同参数对虚警率与检测概率的影响理解似然比阈值与显著性水平的关系并借助蒙特卡洛仿真评估检测性能为雷达、通信等场景下的信号检测问题提供可复用的代码框架与排错思路。1. 从一段“跑不通”的检测代码说起GLRT 与似然比检验到底在解决什么如果你手头有一段信号检测或统计判决的代码文件名里带着GLRT、Statistical test、似然比这些词却不确定它到底在算什么、参数怎么调、结果怎么判那这篇笔记就是写给你的。广义似然比检验GLRT和似然比检验LRT是统计检测理论里最核心的一对工具前者解决“参数未知”的现实问题后者是理论上的最优判决器。很多工程场景——雷达目标检测、通信信号判决、故障诊断、A/B 实验显著性判断——底层都在用这套逻辑。它适合两类人一是刚接触统计检测、想搞懂“为什么不是直接比大小”的新手二是已经在用但被虚警率、门限、未知参数搞得头大的从业者。接下来我会把理论立住再给可复现的 Python 实现最后把踩过的坑摊开讲。2. 似然比检验的判决逻辑从 Neyman-Pearson 到可执行的门限2.1 为什么“比似然比”比“比概率”更靠谱假设你有一串观测数据 (x_1, x_2, \dots, x_N)需要判断它来自 (H_0)只有噪声还是 (H_1)信号加噪声。直觉做法是算两个假设下的概率再比大小但概率密度值本身受量纲和幅度影响直接比会翻车。似然比检验的做法是构造统计量[ \Lambda(\mathbf{x}) \frac{p(\mathbf{x}; H_1)}{p(\mathbf{x}; H_0)} ]然后与门限 (\eta) 比较(\Lambda \eta) 判 (H_1)否则判 (H_0)。Neyman-Pearson 引理告诉我们在给定虚警概率 (P_{FA}) 的前提下这种判决方式能获得最大的检测概率 (P_D)。换句话说它不是“一种”方法而是理论上最优的那一种。工程上我们真正关心的是门限 (\eta) 怎么定才能让虚警率刚好卡在可接受的水平。2.2 门限求解从 (P_{FA}) 反推 (\eta) 的标准步骤门限不是拍脑袋定的而是由虚警率反推。以高斯白噪声中已知信号的检测为例检验统计量在 (H_0) 下服从标准正态分布在 (H_1) 下服从均值为 (d) 的正态分布。给定 (P_{FA} \alpha)门限由正态分布的上分位点决定import numpy as np from scipy.stats import norm def compute_threshold(pfa, N, sigma21.0): 在高斯白噪声、已知信号场景下根据虚警率计算门限。 pfa: 目标虚警概率如 0.01 N: 样本数 sigma2: 噪声方差 返回: 归一化门限 gamma # 检验统计量在 H0 下近似 N(0, sigma2/N) std np.sqrt(sigma2 / N) # 单边检验门限对应 1-pfa 分位点 gamma norm.ppf(1 - pfa) * std return gamma # 示例100 个样本虚警率 1% gamma compute_threshold(0.01, 100) print(f门限 {gamma:.4f})这段代码的逻辑是先确定检验统计量的零假设分布再用分位点函数反推门限。参数pfa直接决定门限高低——pfa越小门限越高虚警少了但检测概率也会降。N越大统计量方差越小门限绝对值越低检测越灵敏。实际使用时sigma2往往未知需要用样本估计这就是 GLRT 要解决的问题。2.3 从 LRT 到 GLRT未知参数怎么处理现实里 (H_1) 下的信号幅度、到达时间、多普勒频率往往未知。这时无法直接写出 (p(\mathbf{x}; H_1))因为它是参数的函数。GLRT 的思路是先用最大似然估计把未知参数估出来再代入似然比[ \Lambda_G(\mathbf{x}) \frac{\max_{\theta_1} p(\mathbf{x}; \theta_1, H_1)}{\max_{\theta_0} p(\mathbf{x}; \theta_0, H_0)} ]如果 (H_0) 下也有未知参数比如噪声方差未知分母同样要做最大似然估计。GLRT 不保证小样本最优但大样本下渐近最优工程上足够用。常见做法是先写出对数似然函数对未知参数求导置零解出估计值再代回构造统计量。下面给一个幅度未知的 GLRT 实现。import numpy as np def glrt_unknown_amplitude(x, s, sigma2): GLRT信号波形 s 已知幅度 A 未知噪声方差 sigma2 已知。 x: 观测向量 s: 信号模板向量 sigma2: 噪声方差 返回: 检验统计量对数似然比 # H1 下 A 的最大似然估计 A_hat np.dot(s, x) / np.dot(s, s) # H1 下最大对数似然 ll1 -np.sum((x - A_hat * s) ** 2) / (2 * sigma2) # H0 下对数似然A0 ll0 -np.sum(x ** 2) / (2 * sigma2) # 对数似然比 lr ll1 - ll0 return lr, A_hat # 模拟信号幅度 0.8噪声标准差 1 np.random.seed(42) N 200 s np.sin(2 * np.pi * 0.1 * np.arange(N)) x 0.8 * s np.random.randn(N) lr, A_hat glrt_unknown_amplitude(x, s, 1.0) print(f对数似然比 {lr:.2f}, 幅度估计 {A_hat:.3f})这里A_hat是匹配滤波器的输出lr可以进一步化简为 (A_hat^2 \cdot |s|^2 / (2\sigma^2))与门限比较即可判决。参数sigma2若未知需要先用 (H_0) 下的样本估计或者联合估计复杂度会上升但思路一致。3. 把 GLRT 跑成可复现的检测流程数据、统计量、判决与性能评估3.1 构造仿真数据让虚警率和检测概率可测理论公式再漂亮不跑数据都是玄学。要验证 GLRT 是否工作必须构造带标签的仿真数据并且能独立控制 (H_0) 和 (H_1) 的生成过程。常见做法是生成纯噪声样本作为 (H_0)生成信号加噪声作为 (H_1)然后分别计算统计量画 ROC 曲线。下面这段代码生成两组数据并计算经验虚警率和检测概率。import numpy as np def simulate_detection(N, num_trials, snr_db, pfa_target): 仿真 GLRT 检测性能。 N: 样本数 num_trials: 每组假设的试验次数 snr_db: 信噪比dB pfa_target: 目标虚警率 返回: 经验虚警率, 经验检测概率 sigma2 1.0 A np.sqrt(2 * sigma2 * 10 ** (snr_db / 10) / N) # 按信噪比反推幅度 s np.sin(2 * np.pi * 0.1 * np.arange(N)) # H0 数据 x0 np.random.randn(num_trials, N) * np.sqrt(sigma2) # H1 数据 x1 A * s np.random.randn(num_trials, N) * np.sqrt(sigma2) # 计算统计量匹配滤波输出平方 stat0 (x0 s) ** 2 / (np.dot(s, s) * sigma2) stat1 (x1 s) ** 2 / (np.dot(s, s) * sigma2) # 根据目标虚警率确定门限 threshold np.percentile(stat0, 100 * (1 - pfa_target)) pfa_emp np.mean(stat0 threshold) pd_emp np.mean(stat1 threshold) return pfa_emp, pd_emp pfa, pd simulate_detection(N100, num_trials10000, snr_db-6, pfa_target0.01) print(f经验虚警率 {pfa:.4f}, 经验检测概率 {pd:.4f})参数说明snr_db控制信号强度-6 dB属于低信噪比场景检测概率通常不高pfa_target是设计目标实际经验虚警率会围绕它波动试验次数越多越接近。threshold用 (H_0) 统计量的分位点确定这是工程上最稳妥的做法——因为实际中 (H_0) 分布往往比理论推导更复杂用数据驱动门限能避免模型失配。3.2 统计量计算中的数值稳定对数域与归一化直接算似然比容易溢出尤其是样本数大时。标准做法是取对数并且把统计量归一化到零假设下均值为 0、方差为 1 的形式。对于高斯场景对数似然比可以化简为匹配滤波输出与门限比较不需要真的算概率密度。但如果噪声分布不是高斯比如拉普拉斯噪声或学生 t 分布就需要老老实实算对数似然。下面是一个通用对数似然比计算框架支持自定义分布。import numpy as np from scipy.stats import norm, laplace def log_likelihood_ratio(x, dist_h0, dist_h1, params_h0, params_h1): 通用对数似然比计算。 x: 观测向量 dist_h0, dist_h1: scipy.stats 分布对象 params_h0, params_h1: 对应分布参数元组 返回: 对数似然比 ll0 np.sum(dist_h0.logpdf(x, *params_h0)) ll1 np.sum(dist_h1.logpdf(x, *params_h1)) return ll1 - ll0 # 高斯噪声 vs 拉普拉斯噪声下的检测 np.random.seed(0) x np.random.randn(50) 0.5 lr_gauss log_likelihood_ratio(x, norm, norm, (0, 1), (0.5, 1)) lr_laplace log_likelihood_ratio(x, laplace, laplace, (0, 1), (0.5, 1)) print(f高斯假设下 LLR {lr_gauss:.2f}) print(f拉普拉斯假设下 LLR {lr_laplace:.2f})注意logpdf在极端值下可能返回-inf实际使用时要加保护比如截断到-1e10。另外如果 (H_0) 和 (H_1) 的噪声方差不同参数元组要相应调整。这个框架的好处是换分布只需换对象和参数不用重写判决逻辑。3.3 性能评估ROC 曲线与检测门限的选取门限选完后必须用 ROC 曲线评估整体性能。ROC 曲线横轴是虚警率纵轴是检测概率曲线越靠左上越好。工程上常看两个指标给定 (P_{FA}) 下的 (P_D)以及达到给定 (P_D) 所需的信噪比。下面代码画 ROC 曲线并计算曲线下面积。import numpy as np import matplotlib.pyplot as plt from sklearn.metrics import roc_curve, auc def plot_roc(stat_h0, stat_h1): 绘制 ROC 曲线并计算 AUC。 stat_h0: H0 下的检验统计量数组 stat_h1: H1 下的检验统计量数组 labels np.concatenate([np.zeros(len(stat_h0)), np.ones(len(stat_h1))]) scores np.concatenate([stat_h0, stat_h1]) fpr, tpr, thresholds roc_curve(labels, scores) roc_auc auc(fpr, tpr) plt.figure() plt.plot(fpr, tpr, labelfAUC {roc_auc:.3f}) plt.plot([0, 1], [0, 1], k--) plt.xlabel(虚警率) plt.ylabel(检测概率) plt.legend() plt.title(GLRT 检测 ROC 曲线) plt.show() return roc_auc # 生成数据并画图 np.random.seed(1) N 100 s np.sin(2 * np.pi * 0.1 * np.arange(N)) x0 np.random.randn(5000, N) x1 0.3 * s np.random.randn(5000, N) stat0 (x0 s) ** 2 / np.dot(s, s) stat1 (x1 s) ** 2 / np.dot(s, s) auc_val plot_roc(stat0, stat1) print(fAUC {auc_val:.3f})roc_curve会自动遍历所有可能门限返回对应的虚警率和检测概率。auc越接近 1 说明检测器区分能力越强。实际调参时如果 AUC 不理想优先检查信号模板是否匹配、噪声方差估计是否准确、样本数是否足够。4. GLRT 落地避坑从虚警率失控到数值溢出的 5 个血泪教训4.1 现象实测虚警率远高于设计值 → 原因噪声方差估计偏小 → 解决用稳健估计或留出纯噪声段设计 (P_{FA}0.01)实测跑到 0.1这是最常见翻车。根因通常是 (H_0) 下的噪声方差被低估导致门限算低了。比如用含信号的样本估计噪声方差自然偏大不恰恰相反如果用了 (H_1) 数据估计噪声方差会偏大门限偏高虚警率偏低。真正导致虚警率飙升的是噪声模型假设为高斯实际有脉冲成分或者样本间相关性被忽略。解决办法是留出一段确认的纯噪声数据专门估计方差或者改用基于排序的稳健估计比如中位数绝对偏差。4.2 现象低信噪比下检测概率骤降 → 原因未知参数估计方差大 → 解决增加样本数或引入先验GLRT 在低信噪比下幅度估计 (A_hat) 的方差很大导致统计量在 (H_1) 下分布展宽与 (H_0) 分布重叠严重。这不是代码 bug是理论极限。能做的增加样本数 (N)因为估计方差按 (1/N) 下降或者如果知道幅度大致范围改用贝叶斯方法引入先验再或者换用更高效的检测器比如能量检测在完全未知时反而更稳。4.3 现象对数似然比算出 NaN 或 inf → 原因概率密度为零或数值下溢 → 解决对数域计算加截断直接算pdf再取对数遇到极端值时pdf返回 0log(0)就是-inf。标准做法是全程用logpdf并且对返回值做截断np.clip(ll, -1e10, 0)。另外如果样本维度很高连logpdf都可能下溢这时要手动展开对数似然表达式把求和拆成逐项计算避免连乘。4.4 现象门限在不同数据集上波动大 → 原因用测试集本身确定门限 → 解决门限必须在独立噪声集上定有人图省事直接在当前数据上取分位点当门限然后报告虚警率。这等于用答案考试虚警率当然准但换一批数据就崩。正确做法门限必须在独立的、确认的 (H_0) 数据集上确定然后固定门限去测新数据。如果无法获得独立噪声集至少要用交叉验证的方式把数据分成定门限和评估两部分。4.5 现象GLRT 统计量与 LRT 差距大 → 原因未知参数过多或样本不足 → 解决检查渐近条件是否满足GLRT 的大样本最优性依赖渐近正态性。如果未知参数个数相对于样本数太多比如 10 个样本估 5 个参数GLRT 性能会明显差于理论 LRT。排查方法固定其他条件逐步增加样本数看 GLRT 性能是否收敛到 LRT。如果不收敛说明模型本身有问题或者参数不可辨识。工程上一般要求样本数至少是未知参数个数的 10 倍以上。5. 进阶技巧用广义似然比做模型阶数选择与变点检测GLRT 不只用于二选一检测还能推广到多假设和变点检测。一个实用技巧是把模型阶数选择转化为似然比检验。比如判断信号是单频还是双频可以构造两个嵌套模型用 GLRT 统计量 (2(\ell_1 - \ell_0)) 与卡方分布分位点比较。下面代码演示如何用 GLRT 做变点检测——判断一段数据中间是否存在均值跳变。import numpy as np from scipy.stats import chi2 def glrt_change_point(x, min_seg10): 单变点检测假设数据在某个位置均值发生跳变。 x: 观测序列 min_seg: 每段最小长度 返回: 最佳变点位置, GLRT 统计量, p 值 N len(x) best_tau, best_stat None, -np.inf for tau in range(min_seg, N - min_seg): seg1, seg2 x[:tau], x[tau:] # 零假设全局均值 mu0 np.mean(x) ll0 -0.5 * np.sum((x - mu0) ** 2) # 备择假设两段各自均值 mu1, mu2 np.mean(seg1), np.mean(seg2) ll1 -0.5 * (np.sum((seg1 - mu1) ** 2) np.sum((seg2 - mu2) ** 2)) stat 2 * (ll1 - ll0) if stat best_stat: best_stat, best_tau stat, tau # 自由度均值从 1 个变 2 个增加 1 p_value 1 - chi2.cdf(best_stat, df1) return best_tau, best_stat, p_value # 模拟前 50 个点均值 0后 50 个点均值 1 np.random.seed(7) x np.concatenate([np.random.randn(50), np.random.randn(50) 1]) tau, stat, p glrt_change_point(x) print(f变点位置 {tau}, GLRT 统计量 {stat:.2f}, p 值 {p:.4f})这段代码遍历所有可能变点对每个位置计算两段均值模型相对于全局均值模型的似然比。best_stat最大处就是最可能的变点。p_value用卡方分布计算自由度等于备择假设比零假设多出的参数个数。注意min_seg防止段太短导致过拟合。实际使用时如果噪声方差未知需要先估计方差再标准化统计量否则卡方近似不准。我自己的习惯是任何 GLRT 代码上线前先用纯噪声跑 10000 次看经验虚警率是否落在设计值附近再用已知信号跑一遍看检测概率是否与理论曲线吻合。这两步过了才敢往真实数据上套。希望帮到你。本文还有配套的精品资源点击获取
广义似然比检验GLRT实战：从理论推导到Python代码实现与避坑指南
简介这份资源是面向统计信号处理学习者与科研人员的MATLAB仿真代码包聚焦广义最大似然比检验GLRT这一假设检验方法帮助理解弱信号与噪声环境下的检测原理。压缩包共7个文件全部为m脚本整体约3KB涵盖数据生成、似然函数计算、临界值确定、决策规则与性能评估等模块可完整复现GLRT的仿真流程。已有452人学习下载适合希望从理论走向实践、掌握检测器设计的中高级读者。通过运行这些脚本读者能直观观察不同参数对虚警率与检测概率的影响理解似然比阈值与显著性水平的关系并借助蒙特卡洛仿真评估检测性能为雷达、通信等场景下的信号检测问题提供可复用的代码框架与排错思路。1. 从一段“跑不通”的检测代码说起GLRT 与似然比检验到底在解决什么如果你手头有一段信号检测或统计判决的代码文件名里带着GLRT、Statistical test、似然比这些词却不确定它到底在算什么、参数怎么调、结果怎么判那这篇笔记就是写给你的。广义似然比检验GLRT和似然比检验LRT是统计检测理论里最核心的一对工具前者解决“参数未知”的现实问题后者是理论上的最优判决器。很多工程场景——雷达目标检测、通信信号判决、故障诊断、A/B 实验显著性判断——底层都在用这套逻辑。它适合两类人一是刚接触统计检测、想搞懂“为什么不是直接比大小”的新手二是已经在用但被虚警率、门限、未知参数搞得头大的从业者。接下来我会把理论立住再给可复现的 Python 实现最后把踩过的坑摊开讲。2. 似然比检验的判决逻辑从 Neyman-Pearson 到可执行的门限2.1 为什么“比似然比”比“比概率”更靠谱假设你有一串观测数据 (x_1, x_2, \dots, x_N)需要判断它来自 (H_0)只有噪声还是 (H_1)信号加噪声。直觉做法是算两个假设下的概率再比大小但概率密度值本身受量纲和幅度影响直接比会翻车。似然比检验的做法是构造统计量[ \Lambda(\mathbf{x}) \frac{p(\mathbf{x}; H_1)}{p(\mathbf{x}; H_0)} ]然后与门限 (\eta) 比较(\Lambda \eta) 判 (H_1)否则判 (H_0)。Neyman-Pearson 引理告诉我们在给定虚警概率 (P_{FA}) 的前提下这种判决方式能获得最大的检测概率 (P_D)。换句话说它不是“一种”方法而是理论上最优的那一种。工程上我们真正关心的是门限 (\eta) 怎么定才能让虚警率刚好卡在可接受的水平。2.2 门限求解从 (P_{FA}) 反推 (\eta) 的标准步骤门限不是拍脑袋定的而是由虚警率反推。以高斯白噪声中已知信号的检测为例检验统计量在 (H_0) 下服从标准正态分布在 (H_1) 下服从均值为 (d) 的正态分布。给定 (P_{FA} \alpha)门限由正态分布的上分位点决定import numpy as np from scipy.stats import norm def compute_threshold(pfa, N, sigma21.0): 在高斯白噪声、已知信号场景下根据虚警率计算门限。 pfa: 目标虚警概率如 0.01 N: 样本数 sigma2: 噪声方差 返回: 归一化门限 gamma # 检验统计量在 H0 下近似 N(0, sigma2/N) std np.sqrt(sigma2 / N) # 单边检验门限对应 1-pfa 分位点 gamma norm.ppf(1 - pfa) * std return gamma # 示例100 个样本虚警率 1% gamma compute_threshold(0.01, 100) print(f门限 {gamma:.4f})这段代码的逻辑是先确定检验统计量的零假设分布再用分位点函数反推门限。参数pfa直接决定门限高低——pfa越小门限越高虚警少了但检测概率也会降。N越大统计量方差越小门限绝对值越低检测越灵敏。实际使用时sigma2往往未知需要用样本估计这就是 GLRT 要解决的问题。2.3 从 LRT 到 GLRT未知参数怎么处理现实里 (H_1) 下的信号幅度、到达时间、多普勒频率往往未知。这时无法直接写出 (p(\mathbf{x}; H_1))因为它是参数的函数。GLRT 的思路是先用最大似然估计把未知参数估出来再代入似然比[ \Lambda_G(\mathbf{x}) \frac{\max_{\theta_1} p(\mathbf{x}; \theta_1, H_1)}{\max_{\theta_0} p(\mathbf{x}; \theta_0, H_0)} ]如果 (H_0) 下也有未知参数比如噪声方差未知分母同样要做最大似然估计。GLRT 不保证小样本最优但大样本下渐近最优工程上足够用。常见做法是先写出对数似然函数对未知参数求导置零解出估计值再代回构造统计量。下面给一个幅度未知的 GLRT 实现。import numpy as np def glrt_unknown_amplitude(x, s, sigma2): GLRT信号波形 s 已知幅度 A 未知噪声方差 sigma2 已知。 x: 观测向量 s: 信号模板向量 sigma2: 噪声方差 返回: 检验统计量对数似然比 # H1 下 A 的最大似然估计 A_hat np.dot(s, x) / np.dot(s, s) # H1 下最大对数似然 ll1 -np.sum((x - A_hat * s) ** 2) / (2 * sigma2) # H0 下对数似然A0 ll0 -np.sum(x ** 2) / (2 * sigma2) # 对数似然比 lr ll1 - ll0 return lr, A_hat # 模拟信号幅度 0.8噪声标准差 1 np.random.seed(42) N 200 s np.sin(2 * np.pi * 0.1 * np.arange(N)) x 0.8 * s np.random.randn(N) lr, A_hat glrt_unknown_amplitude(x, s, 1.0) print(f对数似然比 {lr:.2f}, 幅度估计 {A_hat:.3f})这里A_hat是匹配滤波器的输出lr可以进一步化简为 (A_hat^2 \cdot |s|^2 / (2\sigma^2))与门限比较即可判决。参数sigma2若未知需要先用 (H_0) 下的样本估计或者联合估计复杂度会上升但思路一致。3. 把 GLRT 跑成可复现的检测流程数据、统计量、判决与性能评估3.1 构造仿真数据让虚警率和检测概率可测理论公式再漂亮不跑数据都是玄学。要验证 GLRT 是否工作必须构造带标签的仿真数据并且能独立控制 (H_0) 和 (H_1) 的生成过程。常见做法是生成纯噪声样本作为 (H_0)生成信号加噪声作为 (H_1)然后分别计算统计量画 ROC 曲线。下面这段代码生成两组数据并计算经验虚警率和检测概率。import numpy as np def simulate_detection(N, num_trials, snr_db, pfa_target): 仿真 GLRT 检测性能。 N: 样本数 num_trials: 每组假设的试验次数 snr_db: 信噪比dB pfa_target: 目标虚警率 返回: 经验虚警率, 经验检测概率 sigma2 1.0 A np.sqrt(2 * sigma2 * 10 ** (snr_db / 10) / N) # 按信噪比反推幅度 s np.sin(2 * np.pi * 0.1 * np.arange(N)) # H0 数据 x0 np.random.randn(num_trials, N) * np.sqrt(sigma2) # H1 数据 x1 A * s np.random.randn(num_trials, N) * np.sqrt(sigma2) # 计算统计量匹配滤波输出平方 stat0 (x0 s) ** 2 / (np.dot(s, s) * sigma2) stat1 (x1 s) ** 2 / (np.dot(s, s) * sigma2) # 根据目标虚警率确定门限 threshold np.percentile(stat0, 100 * (1 - pfa_target)) pfa_emp np.mean(stat0 threshold) pd_emp np.mean(stat1 threshold) return pfa_emp, pd_emp pfa, pd simulate_detection(N100, num_trials10000, snr_db-6, pfa_target0.01) print(f经验虚警率 {pfa:.4f}, 经验检测概率 {pd:.4f})参数说明snr_db控制信号强度-6 dB属于低信噪比场景检测概率通常不高pfa_target是设计目标实际经验虚警率会围绕它波动试验次数越多越接近。threshold用 (H_0) 统计量的分位点确定这是工程上最稳妥的做法——因为实际中 (H_0) 分布往往比理论推导更复杂用数据驱动门限能避免模型失配。3.2 统计量计算中的数值稳定对数域与归一化直接算似然比容易溢出尤其是样本数大时。标准做法是取对数并且把统计量归一化到零假设下均值为 0、方差为 1 的形式。对于高斯场景对数似然比可以化简为匹配滤波输出与门限比较不需要真的算概率密度。但如果噪声分布不是高斯比如拉普拉斯噪声或学生 t 分布就需要老老实实算对数似然。下面是一个通用对数似然比计算框架支持自定义分布。import numpy as np from scipy.stats import norm, laplace def log_likelihood_ratio(x, dist_h0, dist_h1, params_h0, params_h1): 通用对数似然比计算。 x: 观测向量 dist_h0, dist_h1: scipy.stats 分布对象 params_h0, params_h1: 对应分布参数元组 返回: 对数似然比 ll0 np.sum(dist_h0.logpdf(x, *params_h0)) ll1 np.sum(dist_h1.logpdf(x, *params_h1)) return ll1 - ll0 # 高斯噪声 vs 拉普拉斯噪声下的检测 np.random.seed(0) x np.random.randn(50) 0.5 lr_gauss log_likelihood_ratio(x, norm, norm, (0, 1), (0.5, 1)) lr_laplace log_likelihood_ratio(x, laplace, laplace, (0, 1), (0.5, 1)) print(f高斯假设下 LLR {lr_gauss:.2f}) print(f拉普拉斯假设下 LLR {lr_laplace:.2f})注意logpdf在极端值下可能返回-inf实际使用时要加保护比如截断到-1e10。另外如果 (H_0) 和 (H_1) 的噪声方差不同参数元组要相应调整。这个框架的好处是换分布只需换对象和参数不用重写判决逻辑。3.3 性能评估ROC 曲线与检测门限的选取门限选完后必须用 ROC 曲线评估整体性能。ROC 曲线横轴是虚警率纵轴是检测概率曲线越靠左上越好。工程上常看两个指标给定 (P_{FA}) 下的 (P_D)以及达到给定 (P_D) 所需的信噪比。下面代码画 ROC 曲线并计算曲线下面积。import numpy as np import matplotlib.pyplot as plt from sklearn.metrics import roc_curve, auc def plot_roc(stat_h0, stat_h1): 绘制 ROC 曲线并计算 AUC。 stat_h0: H0 下的检验统计量数组 stat_h1: H1 下的检验统计量数组 labels np.concatenate([np.zeros(len(stat_h0)), np.ones(len(stat_h1))]) scores np.concatenate([stat_h0, stat_h1]) fpr, tpr, thresholds roc_curve(labels, scores) roc_auc auc(fpr, tpr) plt.figure() plt.plot(fpr, tpr, labelfAUC {roc_auc:.3f}) plt.plot([0, 1], [0, 1], k--) plt.xlabel(虚警率) plt.ylabel(检测概率) plt.legend() plt.title(GLRT 检测 ROC 曲线) plt.show() return roc_auc # 生成数据并画图 np.random.seed(1) N 100 s np.sin(2 * np.pi * 0.1 * np.arange(N)) x0 np.random.randn(5000, N) x1 0.3 * s np.random.randn(5000, N) stat0 (x0 s) ** 2 / np.dot(s, s) stat1 (x1 s) ** 2 / np.dot(s, s) auc_val plot_roc(stat0, stat1) print(fAUC {auc_val:.3f})roc_curve会自动遍历所有可能门限返回对应的虚警率和检测概率。auc越接近 1 说明检测器区分能力越强。实际调参时如果 AUC 不理想优先检查信号模板是否匹配、噪声方差估计是否准确、样本数是否足够。4. GLRT 落地避坑从虚警率失控到数值溢出的 5 个血泪教训4.1 现象实测虚警率远高于设计值 → 原因噪声方差估计偏小 → 解决用稳健估计或留出纯噪声段设计 (P_{FA}0.01)实测跑到 0.1这是最常见翻车。根因通常是 (H_0) 下的噪声方差被低估导致门限算低了。比如用含信号的样本估计噪声方差自然偏大不恰恰相反如果用了 (H_1) 数据估计噪声方差会偏大门限偏高虚警率偏低。真正导致虚警率飙升的是噪声模型假设为高斯实际有脉冲成分或者样本间相关性被忽略。解决办法是留出一段确认的纯噪声数据专门估计方差或者改用基于排序的稳健估计比如中位数绝对偏差。4.2 现象低信噪比下检测概率骤降 → 原因未知参数估计方差大 → 解决增加样本数或引入先验GLRT 在低信噪比下幅度估计 (A_hat) 的方差很大导致统计量在 (H_1) 下分布展宽与 (H_0) 分布重叠严重。这不是代码 bug是理论极限。能做的增加样本数 (N)因为估计方差按 (1/N) 下降或者如果知道幅度大致范围改用贝叶斯方法引入先验再或者换用更高效的检测器比如能量检测在完全未知时反而更稳。4.3 现象对数似然比算出 NaN 或 inf → 原因概率密度为零或数值下溢 → 解决对数域计算加截断直接算pdf再取对数遇到极端值时pdf返回 0log(0)就是-inf。标准做法是全程用logpdf并且对返回值做截断np.clip(ll, -1e10, 0)。另外如果样本维度很高连logpdf都可能下溢这时要手动展开对数似然表达式把求和拆成逐项计算避免连乘。4.4 现象门限在不同数据集上波动大 → 原因用测试集本身确定门限 → 解决门限必须在独立噪声集上定有人图省事直接在当前数据上取分位点当门限然后报告虚警率。这等于用答案考试虚警率当然准但换一批数据就崩。正确做法门限必须在独立的、确认的 (H_0) 数据集上确定然后固定门限去测新数据。如果无法获得独立噪声集至少要用交叉验证的方式把数据分成定门限和评估两部分。4.5 现象GLRT 统计量与 LRT 差距大 → 原因未知参数过多或样本不足 → 解决检查渐近条件是否满足GLRT 的大样本最优性依赖渐近正态性。如果未知参数个数相对于样本数太多比如 10 个样本估 5 个参数GLRT 性能会明显差于理论 LRT。排查方法固定其他条件逐步增加样本数看 GLRT 性能是否收敛到 LRT。如果不收敛说明模型本身有问题或者参数不可辨识。工程上一般要求样本数至少是未知参数个数的 10 倍以上。5. 进阶技巧用广义似然比做模型阶数选择与变点检测GLRT 不只用于二选一检测还能推广到多假设和变点检测。一个实用技巧是把模型阶数选择转化为似然比检验。比如判断信号是单频还是双频可以构造两个嵌套模型用 GLRT 统计量 (2(\ell_1 - \ell_0)) 与卡方分布分位点比较。下面代码演示如何用 GLRT 做变点检测——判断一段数据中间是否存在均值跳变。import numpy as np from scipy.stats import chi2 def glrt_change_point(x, min_seg10): 单变点检测假设数据在某个位置均值发生跳变。 x: 观测序列 min_seg: 每段最小长度 返回: 最佳变点位置, GLRT 统计量, p 值 N len(x) best_tau, best_stat None, -np.inf for tau in range(min_seg, N - min_seg): seg1, seg2 x[:tau], x[tau:] # 零假设全局均值 mu0 np.mean(x) ll0 -0.5 * np.sum((x - mu0) ** 2) # 备择假设两段各自均值 mu1, mu2 np.mean(seg1), np.mean(seg2) ll1 -0.5 * (np.sum((seg1 - mu1) ** 2) np.sum((seg2 - mu2) ** 2)) stat 2 * (ll1 - ll0) if stat best_stat: best_stat, best_tau stat, tau # 自由度均值从 1 个变 2 个增加 1 p_value 1 - chi2.cdf(best_stat, df1) return best_tau, best_stat, p_value # 模拟前 50 个点均值 0后 50 个点均值 1 np.random.seed(7) x np.concatenate([np.random.randn(50), np.random.randn(50) 1]) tau, stat, p glrt_change_point(x) print(f变点位置 {tau}, GLRT 统计量 {stat:.2f}, p 值 {p:.4f})这段代码遍历所有可能变点对每个位置计算两段均值模型相对于全局均值模型的似然比。best_stat最大处就是最可能的变点。p_value用卡方分布计算自由度等于备择假设比零假设多出的参数个数。注意min_seg防止段太短导致过拟合。实际使用时如果噪声方差未知需要先估计方差再标准化统计量否则卡方近似不准。我自己的习惯是任何 GLRT 代码上线前先用纯噪声跑 10000 次看经验虚警率是否落在设计值附近再用已知信号跑一遍看检测概率是否与理论曲线吻合。这两步过了才敢往真实数据上套。希望帮到你。本文还有配套的精品资源点击获取
