Python实战从孟德尔豌豆到A/B测试的似然比检验全解析1. 似然比检验的核心思想似然比检验Likelihood Ratio Test, LRT是统计学中一种强大而通用的假设检验方法它通过比较两种统计模型对数据的解释能力来判断哪个模型更优。这种方法由统计学家奈曼和E.皮尔逊在1928年提出已成为现代统计推断的基石之一。基本原理似然比检验的核心是比较在零假设H₀和备择假设H₁下数据的似然函数最大值之比。这个比值越大说明备择假设比零假设更能解释观察到的数据。数学表达式为 Λ(x₁,x₂,...,xₙ) [sup(θ∈Θ) p(x₁,...,xₙ;θ)] / [sup(θ∈Θ₀) p(x₁,...,xₙ;θ)]其中分子是在全参数空间Θ下的最大似然分母是在零假设参数空间Θ₀下的最大似然在实际应用中我们通常使用对数似然比统计量 -2lnΛ ≈ χ²(df)这个统计量在大样本下服从卡方分布自由度df等于参数空间维度的差异。2. 经典案例孟德尔豌豆实验的Python实现让我们用Python重现孟德尔的经典豌豆实验分析。孟德尔观察了556颗豌豆按颜色和形状分为四类类别观察数理论比例黄圆3159/16绿圆1083/16黄皱1013/16绿皱321/16Python实现步骤import numpy as np from scipy.stats import chi2 # 观察数据 observed np.array([315, 108, 101, 32]) expected_ratio np.array([9, 3, 3, 1]) / 16 expected expected_ratio * observed.sum() # 计算卡方统计量 chi_sq np.sum((observed - expected)**2 / expected) # 计算p值 df len(observed) - 1 p_value 1 - chi2.cdf(chi_sq, df) print(f卡方统计量: {chi_sq:.4f}) print(fP值: {p_value:.4f})输出结果分析卡方统计量0.47P值0.9254这个结果与孟德尔的原始分析一致P值远大于常用的显著性水平0.05因此不能拒绝零假设即观察数据与9:3:3:1的理论比例吻合。3. 现代应用A/B测试中的似然比检验在互联网行业A/B测试是评估产品变更效果的黄金标准。让我们看一个电商网站转化率测试的例子假设我们进行了两种页面设计的测试版本A旧版展示给1000用户150人转化版本B新版展示给1100用户180人转化传统卡方检验 我们可以使用2×2列联表的卡方检验from scipy.stats import chi2_contingency # 构建列联表 table np.array([[150, 850], # A组转化与非转化 [180, 920]]) # B组转化与非转化 chi2_stat, p_val, dof, expected chi2_contingency(table, correctionFalse) print(f卡方统计量: {chi2_stat:.4f}) print(fP值: {p_val:.4f})似然比检验实现 更精确的方法是构建广义线性模型GLM并进行似然比检验import statsmodels.api as sm # 准备数据 n_a, conv_a 1000, 150 n_b, conv_b 1100, 180 # 构建设计矩阵和响应变量 X np.array([0]*n_a [1]*n_b) # 0A组1B组 y np.array([1]*conv_a [0]*(n_a-conv_a) [1]*conv_b [0]*(n_b-conv_b)) # 拟合有约束模型仅截距 model_null sm.GLM(y, np.ones_like(X), familysm.families.Binomial()).fit() # 拟合无约束模型截距分组变量 model_full sm.GLM(y, sm.add_constant(X), familysm.families.Binomial()).fit() # 似然比检验 lr_stat -2*(model_null.llf - model_full.llf) p_value chi2.sf(lr_stat, 1) # 自由度为1 print(f似然比统计量: {lr_stat:.4f}) print(fP值: {p_value:.4f})结果解读 两种方法得到的P值通常接近但似然比检验更灵活可以扩展到更复杂的模型比较场景。4. 实战技巧与常见陷阱4.1 样本量规划在进行似然比检验前合理的样本量规划至关重要。对于A/B测试可以使用功效分析from statsmodels.stats.power import tt_ind_solve_power import math # 假设当前转化率15%期望检测到18%的相对提升即17.7%的绝对转化率 effect_size 2 * math.asin(math.sqrt(0.177)) - 2 * math.asin(math.sqrt(0.15)) # 计算达到80%功效所需的样本量 n_per_group tt_ind_solve_power(effect_sizeeffect_size, alpha0.05, power0.8) print(f每组需要样本量: {math.ceil(n_per_group)})4.2 多重检验问题当同时进行多个假设检验时错误发现率会急剧上升。解决方法包括Bonferroni校正将显著性水平α除以检验次数错误发现率FDR控制from statsmodels.stats.multitest import multipletests p_values [0.01, 0.04, 0.03, 0.21, 0.02] reject, adj_p, _, _ multipletests(p_values, methodbonferroni) print(校正后P值:, adj_p) print(是否拒绝:, reject)4.3 连续变量的似然比检验对于连续变量如检验两组均值差异from scipy.stats import ttest_ind # 生成模拟数据 np.random.seed(42) group_a np.random.normal(5, 1, 50) group_b np.random.normal(5.5, 1, 60) # 传统t检验 t_stat, p_val ttest_ind(group_a, group_b) print(ft检验P值: {p_val:.4f}) # 似然比检验实现 from scipy.optimize import minimize def neg_log_lik(params, data): mu, sigma params n len(data) ll -n/2 * np.log(2*np.pi*sigma**2) - 1/(2*sigma**2)*np.sum((data-mu)**2) return -ll # 拟合有约束模型共同均值 cons ({type: eq, fun: lambda x: x[0]-x[1]}) # μ1μ2 result minimize(neg_log_lik, [5,1], args(np.concatenate([group_a, group_b])), constraintscons) ll_null -result.fun # 拟合无约束模型不同均值 result_a minimize(neg_log_lik, [5,1], argsgroup_a) result_b minimize(neg_log_lik, [5,1], argsgroup_b) ll_full -result_a.fun - result_b.fun # 计算检验统计量 lr_stat -2*(ll_null - ll_full) p_value chi2.sf(lr_stat, 1) print(f似然比统计量: {lr_stat:.4f}) print(fP值: {p_value:.4f})5. 高级应用与扩展5.1 混合效应模型中的似然比检验当数据存在层次结构如用户嵌套于地区时混合效应模型更合适import statsmodels.formula.api as smf # 模拟数据 np.random.seed(123) n_groups 5 group np.repeat(range(n_groups), 20) x np.random.normal(sizen_groups*20) y 2 0.5*x np.random.normal(0, 1, n_groups*20) np.repeat(np.random.normal(0, 0.5, n_groups), 20) data pd.DataFrame({y:y, x:x, group:group}) # 拟合随机截距模型 model_random smf.mixedlm(y ~ x, data, groupsdata[group]).fit() # 拟合固定效应模型相当于忽略分组 model_fixed smf.ols(y ~ x, data).fit() # 似然比检验 lr_stat -2*(model_fixed.llf - model_random.llf) p_value chi2.sf(lr_stat, 1) # 比较1个方差参数 print(f似然比统计量: {lr_stat:.4f}) print(fP值: {p_value:.4f})5.2 模型选择中的信息准则除了假设检验信息准则如AIC、BIC也基于似然原理# 比较多个模型 model1 smf.ols(y ~ x, data).fit() model2 smf.ols(y ~ x np.sin(x), data).fit() print(f模型1 AIC: {model1.aic:.2f}, BIC: {model1.bic:.2f}) print(f模型2 AIC: {model2.aic:.2f}, BIC: {model2.bic:.2f})5.3 贝叶斯视角下的模型比较贝叶斯因子可以看作似然比的贝叶斯版本import pymc3 as pm with pm.Model() as model_null: mu pm.Normal(mu, 0, 10) y_obs pm.Normal(y_obs, mumu, observedy) trace_null pm.sample(1000, tune1000) with pm.Model() as model_alt: mu pm.Normal(mu, 0, 10) sigma pm.HalfNormal(sigma, 10) y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy) trace_alt pm.sample(1000, tune1000) # 计算贝叶斯因子近似 from scipy.stats import bayes_factor bf bayes_factor(trace_null, trace_alt) print(f贝叶斯因子: {bf:.2f})6. 性能优化与最佳实践6.1 大样本下的计算优化当样本量极大时直接计算可能效率低下。可以考虑随机子抽样在线算法分布式计算# 使用minibatch处理大数据 def online_lr_test(data_stream, batch_size1000): ll_null, ll_full 0, 0 for batch in data_stream: # 更新模型参数和似然值 pass return -2*(ll_null - ll_full)6.2 可视化诊断良好的可视化能帮助理解检验结果import matplotlib.pyplot as plt import seaborn as sns # 绘制观察值与期望值对比 categories [黄圆, 绿圆, 黄皱, 绿皱] plt.figure(figsize(10,6)) sns.barplot(xcategories, yobserved, colorblue, alpha0.5, label观察值) sns.barplot(xcategories, yexpected, colorred, alpha0.5, label期望值) plt.legend() plt.title(孟德尔豌豆实验观察值与期望值对比) plt.show()6.3 常见误区与解决方案样本量不足导致检验功效低解决方案预先进行功效分析忽略假设条件如卡方检验要求期望频数≥5解决方案合并类别或使用精确检验p值误解p值不是效应大小解决方案同时报告置信区间多重比较增加假阳性风险解决方案校正p值或控制FDR模型误设错误的分布假设解决方案进行模型诊断检验