总变异 SST = ∑(xᵢ - x̄)²
   ├── 组间变异 SSA = ∑nᵢ(x̄ᵢ - x̄)²  (处理效应, "信号")
   └── 组内变异 SSE = ∑∑(xᵢⱼ - x̄ᵢ)²  (误差, "噪声")

from scipy.stats import f_oneway
import numpy as np

# 3 种广告文案的转化率 (连续化: 假设评分 1-100)
np.random.seed(42)
ad_a = np.random.normal(50, 10, 30)  # A 文案
ad_b = np.random.normal(55, 10, 30)  # B 文案
ad_c = np.random.normal(60, 10, 30)  # C 文案

f_stat, p_val = f_oneway(ad_a, ad_b, ad_c)
print(f"ANOVA: F = {f_stat:.3f}, p = {p_val:.4f}")
# 显著 → 至少两组有差异

import statsmodels.api as sm
from statsmodels.formula.api import ols

# 用 statsmodels 看完整 ANOVA 表
df = pd.DataFrame({
    "score": np.concatenate([ad_a, ad_b, ad_c]),
    "ad": (["A"] * 30) + (["B"] * 30) + (["C"] * 30),
})
model = ols("score ~ C(ad)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))

from scipy.stats import ttest_ind
import itertools

# 3 组, 3 对比较
groups = {"A": ad_a, "B": ad_b, "C": ad_c}
pairs = list(itertools.combinations(groups.keys(), 2))
n_comparisons = len(pairs)
alpha_bonf = 0.05 / n_comparisons

print("事后比较 (Bonferroni 修正):")
for g1, g2 in pairs:
    t_stat, p = ttest_ind(groups[g1], groups[g2])
    sig = "显著" if p < alpha_bonf else "不显著"
    print(f"  {g1} vs {g2}: t = {t_stat:.3f}, p = {p:.4f} → {sig}")
# 修正后阈值 0.05/3 = 0.0167

from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(df["score"], df["ad"], alpha=0.05)
print(tukey)
# 输出: 每对比较, 是否拒绝 + 置信区间

# 偏 eta² (控制其他因素)
from statsmodels.stats.anova import anova_lm
result = anova_lm(model, typ=2)
eta_sq = result["sum_sq"]["C(ad)"] / result["sum_sq"].sum()
print(f"η² = {eta_sq:.3f}")

from scipy.stats import levene, bartlett

# Levene 检验 (中位数, 稳健)
stat, p = levene(ad_a, ad_b, ad_c)
print(f"Levene 检验: stat = {stat:.3f}, p = {p:.4f}")
# p > 0.05 → 方差齐性, ANOVA 适用

# Bartlett 检验 (正态假设)
stat, p = bartlett(ad_a, ad_b, ad_c)
print(f"Bartlett 检验: stat = {stat:.3f}, p = {p:.4f}")

# 例: 文案 (3 种) × 时段 (2 种: 早/晚)
import pandas as pd
import numpy as np
np.random.seed(42)
df = pd.DataFrame({
    "ad": np.repeat(["A", "B", "C"], 20),
    "time": (["morning"] * 10 + ["evening"] * 10) * 3,
    "score": np.concatenate([
        np.random.normal(50, 5, 10), np.random.normal(55, 5, 10),  # A
        np.random.normal(55, 5, 10), np.random.normal(60, 5, 10),  # B
        np.random.normal(60, 5, 10), np.random.normal(50, 5, 10),  # C (晚上低)
    ]),
})

# 双因素 ANOVA (含交互)
model = ols("score ~ C(ad) * C(time)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))
# 输出:
#                 sum_sq   df          F    PR(>F)
# C(ad)            ...    2         ...      ...
# C(time)          ...    1         ...      ...
# C(ad):C(time)    ...    2         ...      ...   ← 交互效应
# Residual         ...   54

# 反例: 5 组, 跑 10 次 t 检验, 至少 1 次假阳性的概率
p_at_least_one = 1 - (1 - 0.05)**10
print(f"5 组, 跑 10 次 t 检验, 至少 1 次假阳性概率 = {p_at_least_one:.1%}")
# 40%!

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 模拟 4 种营销活动, 各 50 个用户, 转化评分 (1-100)
np.random.seed(42)
data = pd.DataFrame({
    "campaign": np.repeat(["邮件", "短信", "推送", "广告"], 50),
    "score": np.concatenate([
        np.random.normal(60, 10, 50),
        np.random.normal(65, 10, 50),
        np.random.normal(70, 10, 50),
        np.random.normal(68, 10, 50),
    ]),
})

# 1. 描述统计
print(data.groupby("campaign")["score"].agg(["mean", "std", "count"]))

# 2. 方差齐性
from scipy.stats import levene
groups = [data[data["campaign"] == c]["score"] for c in data["campaign"].unique()]
stat, p_lev = levene(*groups)
print(f"\nLevene 方差齐性: p = {p_lev:.3f}")

# 3. 单因素 ANOVA
model = ols("score ~ C(campaign)", data=data).fit()
print("\nANOVA 表:")
print(sm.stats.anova_lm(model, typ=2))

# 4. 事后检验
tukey = pairwise_tukeyhsd(data["score"], data["campaign"], alpha=0.05)
print("\nTukey HSD 事后检验:")
print(tukey)

# 5. 效应大小
ss_campaign = sm.stats.anova_lm(model, typ=2)["sum_sq"]["C(campaign)"]
ss_total = data["score"].var(ddof=1) * (len(data) - 1)
eta_sq = ss_campaign / ss_total
print(f"\nη² = {eta_sq:.3f}")

# 特征选择: 类别型特征对连续 target 的影响
import pandas as pd
from sklearn.datasets import fetch_california_housing

data = fetch_california_housing()
df = pd.DataFrame(data.data, columns=data.feature_names)
df["target"] = data.target
df["high_income"] = (df["Population"] > df["Population"].median()).astype(str)

# 看"高人口区" vs "低人口区" 房价均值
model = ols("target ~ C(high_income)", data=df.sample(1000, random_state=0)).fit()
print(sm.stats.anova_lm(model, typ=2))
# 显著 → 人口密度跟房价相关