1. 这不是统计学课件，而是一份能直接跑通的Python统计实战手册

“All Statistical Tests using Python: Mastering Statistics | Part — 1”——看到这个标题，我第一反应不是点开看，而是顺手关掉浏览器标签页。原因很简单：过去三年里，我亲手调试过17个声称“全覆盖”的统计Python教程，其中14个卡在 t检验p值计算结果和SPSS不一致 ，2个在 卡方检验的期望频数校验逻辑上出错 ，还有1个干脆把 单样本Wilcoxon符号秩检验误写成配对样本t检验的替代方案 ，导致用户在非正态小样本场景下得出完全错误的结论。这不是理论疏漏，是实操断层。所以今天这篇，不讲中心极限定理的哲学意义，不画正态分布曲线的数学美感，只做一件事： 把每一种基础统计检验，变成你明天就能粘贴进Jupyter Notebook、输入真实数据、立刻得到可解释结果的可靠代码块 。核心关键词—— Python统计检验、scipy.stats实操、假设检验落地、p值解读陷阱、小样本校验逻辑 ——全部嵌入在真实调试场景中。适合三类人：刚转行的数据分析师（需要快速交付可信报告）、生物/心理/教育等领域的科研人员（常被审稿人追问检验前提是否满足）、以及被“统计显著”四个字反复折磨的产品经理（终于能看懂AB测试背后的数字到底在说什么）。它不承诺让你成为统计学家，但能确保你下次提交分析结果时，心里有底。

2. 整体设计思路：为什么放弃“教科书式分类”，选择“问题驱动流程图”

2.1 拒绝按检验名称罗列，转向“你的数据长什么样？”的决策树

市面上90%的统计教程，结构都是“第1章 t检验 → 第2章 卡方检验 → 第3章 ANOVA……”。这种结构在考试中很高效，但在真实项目里是灾难。上周帮一个电商团队分析促销效果，他们原始数据是：3个不同折扣力度（5折、7折、8.5折）下，各15天的订单转化率（百分比）。按传统目录，你会先翻到“ANOVA”章节——但问题来了：转化率是比例数据，严格来说不满足ANOVA对连续正态分布的假设；而且3组样本量都只有15，小样本下正态性检验本身就不稳定。这时候，教科书式的路径会把你带进死胡同：要么强行用ANOVA被质疑方法论，要么卡在“该用什么非参数检验”上反复查资料。我们的设计反其道而行： 所有检验不按名称分组，而按三个硬性条件分流 ：

1. 你的因变量（Y）是什么类型？ 是连续数值（如销售额、响应时间）、分类标签（如“购买/未购买”）、还是有序等级（如满意度1-5分）？
2. 你的自变量（X）有几个水平？ 是两组（A/B）、三组及以上（A/B/C）、还是连续变量（如广告投放金额）？
3. 你的数据结构是什么？ 是独立样本（两组人互不重叠）、配对样本（同一组人测两次）、还是重复测量（同一组人测多次）？

这三个问题的答案，直接对应到一张极简决策表。比如上面电商案例：Y=转化率（连续数值）、X=3个折扣水平（三组）、数据结构=独立样本（每天数据独立）→ 立刻锁定“单因素方差分析（One-way ANOVA）或Kruskal-Wallis H检验”。接下来才进入具体检验的选择逻辑，而不是在10种检验名称里大海捞针。

2.2 工具链锁定scipy.stats：为什么不用statsmodels或pingouin？

选型不是技术炫技，而是为稳定性让路。我对比过scipy.stats、statsmodels、pingouin在2023年Q4的维护状态：

• scipy.stats ：底层C实现，计算速度最快；API极度稳定， scipy.stats.ttest_ind() 从1.2.0版本（2019年）到当前1.11.4（2023年），参数名、返回值结构零变更；文档中每个函数都明确标注“适用前提”和“注意事项”，比如 scipy.stats.chi2_contingency() 会直接警告“当期望频数<5的单元格超过20%，结果可能不可靠”。

• statsmodels ：功能强大，但API迭代激进。2022年 statsmodels.stats.power 模块重构成 statsmodels.stats.api ，导致大量旧脚本报错；其 anova_lm() 函数默认使用Type II SS，而多数用户实际需要Type III，需额外参数配置，新手极易踩坑。

• pingouin ：语法极其简洁（如 pg.ttest() 一行搞定），但底层仍调用scipy，属于封装层。问题在于：它把“假设检验前提验证”全隐藏了。比如 pg.ttest() 不强制要求你先做正态性检验，用户直接喂入偏态数据，函数照样返回p值——这等于把风险转嫁给使用者。

因此，本系列所有代码， 只依赖scipy.stats（+ numpy/pandas做数据预处理） 。不是因为它最炫，而是因为它的错误提示最直白、前提检查最透明、结果复现性最高。当你在生产环境跑一个关键AB测试时，你不需要“优雅”，你需要“确定性”。

2.3 “Part-1”的边界划定：聚焦高频刚需，砍掉学术冷门

标题里“Part-1”不是营销话术，是严格的范围控制。我们只覆盖 工作中出现频率TOP 8、且最容易用错的检验 ，按实际使用强度排序：

1. 单样本t检验 （检验均值是否等于某个理论值，如“用户平均停留时长是否达到行业基准8分钟”）
2. 独立样本t检验 （两组独立人群比较，如“A/B版注册流程的转化率差异”）
3. 配对样本t检验 （同一组人在两个时间点/条件下的比较，如“用户升级APP前后每日使用时长”）
4. 单因素方差分析（ANOVA） （三组及以上独立样本均值比较，如“三种客服话术的客户满意度评分”）
5. Kruskal-Wallis H检验 （ANOVA的非参数替代，当ANOVA前提不满足时启用）
6. 卡方检验（独立性） （两个分类变量是否相关，如“性别与购买品类偏好是否有关联”）
7. Fisher精确检验 （卡方检验的小样本救急方案，当期望频数过低时）
8. Pearson相关系数 （两个连续变量线性相关强度，如“广告曝光量与销售额”）

为什么砍掉曼-惠特尼U检验？因为它和独立样本t检验解决的是同一类问题（两组独立样本），而Kruskal-Wallis是其自然延伸，掌握后者自动覆盖前者。为什么没有多元回归？那是Part-2的内容——Part-1的使命，就是让你把最常遇到的“两组比、三组比、分类关联、相关性”这四类问题，用最稳妥的方式闭环。

3. 核心细节解析：每一个参数背后，都是血泪教训换来的经验

3.1 单样本t检验：你以为的“简单”，藏着最大的p值陷阱

单样本t检验的公式看似简单：t = (x̄ - μ₀) / (s/√n)。但实操中， 90%的错误发生在μ₀的设定和p值解读上 。

先看一个真实案例：某SaaS公司想验证“新用户7日留存率是否达到产品目标值35%”。运营同学导出数据，计算得样本均值x̄=32.1%，标准差s=8.5%，n=120。他直接运行：

from scipy import stats
import numpy as np
# 假设data是120个用户的留存率（小数形式，如0.321）
t_stat, p_value = stats.ttest_1samp(data, popmean=0.35)
print(f"t={t_stat:.3f}, p={p_value:.4f}")

输出：t=-3.921, p=0.0001。他兴奋地写报告：“显著低于目标！”——但问题来了： p=0.0001代表什么？是“低于目标”的概率？还是“等于目标”的概率？

提示：p值永远不是“原假设为真”的概率，而是“在原假设为真时，观察到当前样本或更极端样本的概率”。这里原假设H₀是“真实均值μ=0.35”，所以p=0.0001意味着：如果真实留存率真是35%，那么抽到均值≤32.1%（或≥37.9%，因为双尾）的样本，概率只有0.01%。这是一个极小概率事件，因此我们拒绝H₀，接受备择假设H₁：“真实均值≠35%”。但注意！ 拒绝H₀不等于证明H₁为真，更不等于知道方向 。要确认是“低于”还是“高于”，必须看t统计量的符号：t为负，说明样本均值小于μ₀，因此结论是“显著低于35%”。

更隐蔽的陷阱在 popmean 参数。很多同学把目标值35%直接写成35，而非0.35。 scipy.stats.ttest_1samp() 不会报错，但会返回完全错误的t值（因为数据是0.321，而你传入35，差了100倍）。 我的实操心得：所有比例、百分比数据，在输入前必须统一转换为[0,1]区间，并用注释明确标出 ：

# data: array of 120 users' 7-day retention rates, e.g., [0.321, 0.285, ...]
# popmean=0.35 means "target is 35%", NOT 35
t_stat, p_value = stats.ttest_1samp(data, popmean=0.35, alternative='less') # 关键！指定单尾检验

这里新增了 alternative='less' 。因为业务问题明确是“是否低于目标”，用单尾检验更敏感（p值约为双尾的一半），且避免了“双尾p=0.0001，但实际关心的是单侧”的逻辑断裂。scipy默认双尾，这是必须手动修正的点。

3.2 独立样本t检验：方差齐性不是可选项，而是必经关卡

两组独立样本比较（如A/B测试）， scipy.stats.ttest_ind() 有两个关键参数： equal_var （是否假设方差齐性）和 nan_policy （如何处理缺失值）。 85%的线上事故源于 equal_var 的误设 。

场景还原：某APP做按钮颜色AB测试，A组（蓝色按钮）5000用户，转化率均值12.3%，标准差2.1%；B组（绿色按钮）4800用户，均值13.1%，标准差3.8%。同学直接运行：

t_stat, p_value = stats.ttest_ind(group_a, group_b) # 默认equal_var=True

结果p=0.023，宣布B组显著胜出。但Levene检验（ scipy.stats.levene() ）显示两组方差差异极显著（p<0.001）。这意味着 equal_var=True 的假设被严重违反，此时t检验的I类错误率（假阳性）会飙升——你看到的p=0.023，实际可能是0.08甚至更高。

注意： scipy.stats.ttest_ind() 的 equal_var=True 对应的是 Student's t-test ，要求两总体方差相等； equal_var=False 对应 Welch's t-test ，不假设方差齐性，且对小样本更稳健。现代统计实践已普遍推荐Welch检验作为默认选项，因为方差齐性在现实中极少严格成立。

正确流程必须是三步：

1. 先做方差齐性检验 （Levene或Bartlett）；
2. 根据结果决定t检验类型 ；
3. 明确报告所用检验名称 （避免审稿人质疑）。

from scipy import stats
# Step 1: Test for equal variance
levene_stat, levene_p = stats.levene(group_a, group_b)
print(f"Levene test: statistic={levene_stat:.3f}, p={levene_p:.4f}")

# Step 2: Choose t-test based on result
if levene_p > 0.05:  # 方差齐性成立
    t_stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=True)
    test_name = "Student's t-test"
else:  # 方差不齐，用Welch检验
    t_stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=False)
    test_name = "Welch's t-test"

print(f"Using {test_name}: t={t_stat:.3f}, p={p_value:.4f}")

这个三步法，是我给所有数据分析新人的第一条硬性规范。它多花3行代码，但能避免90%的统计误读。

3.3 配对样本t检验：数据结构错位，是比计算错误更致命的错误

配对检验的核心是“同一主体在两个条件下的测量”，例如：20名用户分别体验旧版和新版搜索框，记录每次的搜索耗时。数据结构必须是 20行×2列 （每行是一个用户，两列是旧版/新版耗时）。但常见错误是把数据整理成 40行×1列 （20个旧版数据+20个新版数据），然后错误地调用 stats.ttest_ind() 。

后果是什么？ ttest_ind() 会把这40个数据当成40个独立样本，完全无视“用户ID”的配对关系，统计效力暴跌，p值失真。 配对检验的威力，恰恰来自消除个体差异 。比如用户A天生手速快（旧版2s，新版1.8s），用户B手速慢（旧版8s，新版7.5s），配对检验关注的是差值（-0.2s, -0.5s），而独立检验却把它们混在40个随机点中。

正确做法： 用 stats.ttest_rel() ，且数据必须是两个等长数组，索引一一对应 。

# 正确：group_old[i] 和 group_new[i] 是同一个用户的两次测量
t_stat, p_value = stats.ttest_rel(group_old, group_new)

# 错误：把所有旧版数据拼接在前，所有新版在后
all_data = np.concatenate([group_old, group_new]) # 绝对禁止！
# ttest_ind(all_data[:20], all_data[20:]) # 这是伪配对，结果无效

我的实操心得：在写代码前，先用 pandas.DataFrame 把数据框出来，肉眼确认前5行：

user_id	old_time	new_time
U001	2.1	1.9
U002	5.3	4.8
...	...	...
只要这一列对齐了， ttest_rel() 就不可能错。这是比背公式更重要的习惯。

3.4 ANOVA与Kruskal-Wallis：何时该“降级”用非参数检验？

单因素ANOVA要求三大前提： 独立性、正态性、方差齐性 。但现实数据往往只满足第一条。比如分析三家门店的月度客单价，样本量各30，但直方图显示明显右偏（大量低价订单，少数高价订单拉高均值）。

此时，强行用 scipy.stats.f_oneway() 会怎样？ANOVA对正态性有一定鲁棒性，但当偏度>2或峰度>4时，I类错误率会从5%升至15%以上。更危险的是， 很多人用ANOVA后不做事后检验（post-hoc），直接说“三组有差异”，却不告诉老板哪两组不同 。

正确路径是“双轨制”：

• 轨道A（参数法） ：先用 scipy.stats.shapiro() 检验每组正态性（注意：shapiro对小样本敏感，n<50才推荐），再用 levene() 检验方差齐性。三者全通过，才用ANOVA + Tukey HSD事后检验。
• 轨道B（非参数法） ：任一前提失败，立即切换Kruskal-Wallis（ scipy.stats.kruskal() ），它只假设数据是独立的、来自连续分布，对形状无要求。

但Kruskal-Wallis也有坑：它检验的是“中位数是否全相等”，而非“均值”。当分布严重偏斜时，中位数和均值差异很大，结论可能和业务直觉冲突。这时， 必须补充效应量（effect size） ，比如 epsilon_squared （Kruskal-Wallis的效应量，范围0-1，>0.2为中等效应）。p值告诉你“是否不同”，效应量告诉你“不同有多大”。

from scipy import stats
import numpy as np

# 假设three_stores是三个数组，各30个客单价
# Step 1: Check normality per group
normal_pvals = [stats.shapiro(group)[1] for group in three_stores]
print("Shapiro p-values:", [f"{p:.4f}" for p in normal_pvals])

# Step 2: If any p < 0.05, skip ANOVA, go to Kruskal-Wallis
h_stat, p_value = stats.kruskal(*three_stores)
print(f"Kruskal-Wallis: H={h_stat:.3f}, p={p_value:.4f}")

# Step 3: Calculate effect size (epsilon-squared)
N = sum(len(g) for g in three_stores)  # total sample size
k = len(three_stores)  # number of groups
epsilon_sq = (h_stat - k + 1) / (N - k)
print(f"Epsilon-squared = {epsilon_sq:.3f} (small<0.01, medium<0.06, large>0.14)")

记住： 没有“更高级”的检验，只有“更匹配数据”的检验 。当数据不服从正态，坚持用ANOVA不是严谨，是固执。

3.5 卡方检验与Fisher精确检验：小样本的生存法则

卡方检验（ scipy.stats.chi2_contingency() ）的致命弱点： 要求每个单元格的期望频数≥5 。但业务数据常是稀疏的。例如分析“用户是否流失”（是/否）与“是否使用过VIP功能”（是/否）的关系，交叉表可能是：

	使用VIP	未使用VIP	总计
流失	3	42	45
未流失	1	158	159
总计	4	200	204

计算期望频数：(45×4)/204 ≈ 0.88 < 5。此时卡方检验的p值不可信。 chi2_contingency() 会返回一个警告，但很多同学直接忽略。

解决方案不是“硬着头皮用”，而是 无缝切换到Fisher精确检验（ scipy.stats.fisher_exact() ） 。它不依赖大样本近似，直接计算所有可能表格的精确概率，是小样本的黄金标准。

import numpy as np
from scipy import stats

# 构建2x2列联表（必须是整数！）
observed = np.array([[3, 42],
                     [1, 158]])

# Step 1: Run chi-square (for reference, but check warning)
chi2, p_chi2, dof, expected = stats.chi2_contingency(observed)
print(f"Chi-square p-value: {p_chi2:.4f}")
print(f"Expected frequencies:\n{expected}")

# Step 2: If any expected < 5, use Fisher
if np.any(expected < 5):
    oddsratio, p_fisher = stats.fisher_exact(observed, alternative='two-sided')
    print(f"Fisher exact p-value: {p_fisher:.4f}")
    print(f"Odds ratio: {oddsratio:.3f}")

注意： fisher_exact() 只支持2x2表。如果是2x3或3x3，需用 scipy.stats.boschloo_exact() （较新函数，需scipy>=1.7.0）或蒙特卡洛模拟。但2x2覆盖了80%的业务场景。

4. 实操全流程：从原始数据到可交付报告的7个硬核步骤

4.1 步骤1：数据加载与探索性分析（EDA）——别跳过这10分钟

所有统计检验的根基，是清晰理解你的数据。我强制要求自己执行以下EDA清单，少一项都不开始建模：

1. 数据概览 ： df.info() 看缺失值、数据类型； df.describe() 看数值型变量的均值、标准差、分位数；
2. 分布可视化 ：对Y变量（如转化率）画直方图+核密度估计（ sns.histplot(data=df, x='conversion_rate', kde=True) ），肉眼判断偏斜、多峰；
3. 分组对比 ：对X变量（如实验组）画箱线图（ sns.boxplot(data=df, x='group', y='conversion_rate') ），看中位数、离散度、异常值；
4. 缺失值处理 ：明确缺失机制（随机缺失？系统缺失？）。若缺失率<5%，用均值/中位数填充；>5%，必须分析缺失是否与Y相关（如高价值用户更不愿填问卷），否则填充会引入偏差；
5. 异常值识别 ：用IQR法（Q1-1.5×IQR, Q3+1.5×IQR）标记， 但绝不直接删除 ！先查证：是录入错误（如转化率120%）？还是真实极端值（如CEO下单1000万）？前者修正，后者保留并注明。

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# 加载数据（示例：ab_test.csv）
df = pd.read_csv('ab_test.csv')
print("Data shape:", df.shape)
print("\nMissing values:")
print(df.isnull().sum())

# EDA for conversion_rate
plt.figure(figsize=(12, 8))

# Histogram + KDE
plt.subplot(2, 2, 1)
sns.histplot(data=df, x='conversion_rate', kde=True, bins=30)
plt.title('Distribution of Conversion Rate')

# Boxplot by group
plt.subplot(2, 2, 2)
sns.boxplot(data=df, x='group', y='conversion_rate')
plt.title('Conversion Rate by Group')

# Q-Q plot for normality check
from scipy import stats
plt.subplot(2, 2, 3)
stats.probplot(df['conversion_rate'], dist="norm", plot=plt)
plt.title('Q-Q Plot for Normality')

# Correlation heatmap (if multiple X)
if 'ad_spend' in df.columns:
    plt.subplot(2, 2, 4)
    sns.heatmap(df[['conversion_rate', 'ad_spend', 'page_views']].corr(), 
                annot=True, cmap='coolwarm')
    plt.title('Correlation Matrix')

plt.tight_layout()
plt.show()

这段代码生成4张图，10分钟内完成数据“体检”。它不产生p值，但决定了你该用t检验还是Mann-Whitney，该用Pearson还是Spearman。 跳过EDA的统计分析，就像没试飞就让飞机上天 。

4.2 步骤2：前提条件验证——写进报告的“方法论声明”

统计检验不是黑箱，你的报告必须向读者（老板、研发、法务）透明展示“为什么选这个检验”。这需要一份标准化的“前提验证表”。

以独立样本t检验为例，报告中必须包含：

检验前提	验证方法	本数据结果	是否满足	依据说明
独立性	数据采集描述	A/B组用户随机分配，无交叉	是	实验设计文档v2.1
正态性（每组）	Shapiro-Wilk检验	Group A: p=0.12; Group B: p=0.08	是	p>0.05，不拒绝正态假设
方差齐性	Levene检验	p=0.21	是	p>0.05，方差齐性成立
样本量（每组）	计数	n_A=5000, n_B=4800	是	n>30，中心极限定理适用

这张表，必须用代码自动生成，而非手动填写。这样既杜绝笔误，又保证可复现。核心是把 shapiro() 、 levene() 的结果，用 pd.DataFrame 结构化输出。

def check_ttest_assumptions(group_a, group_b, alpha=0.05):
    """Automatically generate the assumption check table"""
    from scipy import stats
    import pandas as pd
    
    # Shapiro for each group
    shapiro_a = stats.shapiro(group_a)
    shapiro_b = stats.shapiro(group_b)
    
    # Levene for variance
    levene_stat, levene_p = stats.levene(group_a, group_b)
    
    # Compile results
    results = {
        'Assumption': ['Independence', 'Normality (Group A)', 'Normality (Group B)', 
                       'Equal Variance', 'Sample Size (A)', 'Sample Size (B)'],
        'Method': ['Design description', 'Shapiro-Wilk', 'Shapiro-Wilk', 
                   'Levene test', 'Count', 'Count'],
        'Result': ['N/A', f'p={shapiro_a[1]:.3f}', f'p={shapiro_b[1]:.3f}',
                   f'p={levene_p:.3f}', f'n={len(group_a)}', f'n={len(group_b)}'],
        'Met?': ['Yes', 'Yes' if shapiro_a[1] > alpha else 'No',
                 'Yes' if shapiro_b[1] > alpha else 'No',
                 'Yes' if levene_p > alpha else 'No',
                 'Yes' if len(group_a) >= 30 else 'No',
                 'Yes' if len(group_b) >= 30 else 'No'],
        'Justification': ['Random assignment', 
                          f'Shapiro p>{alpha}, no evidence against normality',
                          f'Shapiro p>{alpha}, no evidence against normality',
                          f'Levene p>{alpha}, no evidence against equal variance',
                          'n>=30 satisfies CLT', 'n>=30 satisfies CLT']
    }
    
    return pd.DataFrame(results)

# Usage
assump_df = check_ttest_assumptions(df[df['group']=='A']['conversion_rate'], 
                                   df[df['group']=='B']['conversion_rate'])
print(assump_df.to_string(index=False))

这份输出，可直接复制进分析报告。它让统计过程从“黑箱操作”变为“透明工程”。

4.3 步骤3：核心检验执行与结果提取——一行代码，三行解读

执行检验只是开始，关键是如何把 scipy.stats 返回的元组，转化为业务语言。以ANOVA为例， f_oneway() 返回 (F_stat, p_value) ，但这远远不够。

你需要：

• F统计量的临界值对比 ：查F分布表，或用 scipy.stats.f.ppf(0.95, dfn, dfd) 计算临界值，说明“F=5.21 > F_critical=3.05，因此拒绝H₀”；
• 效应量η²（Eta-squared） ：衡量X对Y变异的解释比例，公式η² = SS_between / SS_total。 scipy 不直接提供，需手动计算；
• 事后检验（Tukey HSD） ：ANOVA只告诉你“至少有两组不同”，Tukey告诉你“哪两组不同，差异多大”。

from scipy import stats
import numpy as np
from statsmodels.stats.multicomp import pairwise_tukeyhsd

def anova_with_posthoc(groups, alpha=0.05):
    """
    Run ANOVA + Tukey HSD, return comprehensive results
    groups: list of arrays, e.g., [group_a, group_b, group_c]
    """
    # Step 1: ANOVA
    f_stat, p_anova = stats.f_oneway(*groups)
    
    # Step 2: Calculate eta-squared
    # Manual SS calculation
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)
    ss_total = np.sum((all_data - grand_mean) ** 2)
    
    ss_between = 0
    for i, group in enumerate(groups):
        group_mean = np.mean(group)
        ss_between += len(group) * (group_mean - grand_mean) ** 2
    eta_sq = ss_between / ss_total
    
    # Step 3: Tukey HSD
    # Stack data for Tukey (requires 1D array + labels)
    stacked_data = np.concatenate(groups)
    labels = np.concatenate([np.full(len(g), f'Group_{i+1}') for i, g in enumerate(groups)])
    tukey_result = pairwise_tukeyhsd(stacked_data, labels, alpha=alpha)
    
    # Print comprehensive report
    print(f"=== ANOVA Results ===")
    print(f"F-statistic: {f_stat:.4f}")
    print(f"p-value: {p_anova:.4f}")
    print(f"Eta-squared (η²): {eta_sq:.4f} (Small<0.01, Medium<0.06, Large>0.14)")
    print(f"\n=== Tukey HSD Post-hoc ===")
    print(tukey_result)
    
    return {'f_stat': f_stat, 'p_anova': p_anova, 'eta_sq': eta_sq, 'tukey': tukey_result}

# Usage
groups = [df[df['store']=='A']['revenue'], 
          df[df['store']=='B']['revenue'], 
          df[df['store']=='C']['revenue']]
results = anova_with_posthoc(groups)

输出中， Tukey HSD 表格会明确列出：

• Group_1-Group_2 : diff=-1250.3, lower=-1890.1, upper=-610.5, reject=True → A店比B店平均少赚1250元，95%置信区间不包含0，差异显著。

这才是业务方能看懂的语言。

4.4 步骤4：结果可视化——让p值自己说话

数字再精准，不如一张图直观。但可视化不是画个柱状图就完事。针对不同检验，有专属图表：

• t检验/ANOVA ：分组箱线图 + 显著性标记（***表示p<0.001）；
• 相关性 ：散点图 + 回归线 + Pearson r值标注；
• 卡方/Fisher ：堆叠柱状图（显示比例）+ 卡方统计量标注。

def plot_ttest_results(group_a, group_b, group_names=['A', 'B'], 
                      title='Conversion Rate Comparison'):
    """Plot t-test results with significance stars"""
    import matplotlib.pyplot as plt
    import seaborn as sns
    
    # Prepare data for boxplot
    data_long = pd.DataFrame({
        'value': np.concatenate([group_a, group_b]),
        'group': [group_names[0]]*len(group_a) + [group_names[1]]*len(group_b)
    })
    
    plt.figure(figsize=(8, 6))
    ax = sns.boxplot(data=data_long, x='group', y='value')
    
    # Add significance annotation
    t_stat, p_val = stats.ttest_ind(group_a, group_b, equal_var=False)
    if p_val < 0.001:
        sig_star = '***'
    elif p_val < 0.01:
        sig_star = '**'
    elif p_val < 0.05:
        sig_star = '*'
    else:
        sig_star = 'ns'
    
    # Draw line and star between groups
    y_max = data_long['value'].max() * 1.05
    plt.plot([0, 0, 1, 1], [y_max, y_max*1.02, y_max*1.02, y_max], 
             lw=1.5, c='black')
    plt.text(0.5, y_max*1.03, sig_star, ha='center', va='bottom', fontsize=14)
    
    plt.title(f'{title}\n(p={p_val:.4f})', fontsize=14)
    plt.ylabel('Conversion Rate (%)')
    plt.show()

# Usage
plot_ttest_results(df[df['group']=='A']['conversion_rate'], 
                   df[df['group']=='B']['conversion_rate'])

这张图，老板扫一眼就知道：“B组箱子整体更高，顶上有个***，说明差异非常显著”。可视化是统计分析的最终翻译器。

4.5 步骤5：报告撰写模板——把技术语言翻译成业务决策

分析做完，报告是交付物。我用固定模板，确保信息无遗漏：

1. 背景与问题 （1句话）：“为评估新版注册流程对转化率的影响，对比A/B两组7日数据。”
2. 数据概况 （3行）：“共收集有效用户12,500人（A组6,300；B组6,200）。转化率均值A=11.2%±2.3%，B=12.8%±2.1%。”
3. 方法论 （2行）：“采用Welch's t-test（因Levene检验p=0.03，方差不齐）。检验水准α=0.05，双尾。”
4. 核心结果 （2行）：“B