[数据分析]A/A测试：为什么它是A/B测试前的必要步骤？

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import statsmodels.api as sm
from statsmodels.stats import power, proportion

# 设置随机种子确保结果可重现
np.random.seed(42)

# 生成模拟用户数据
n_users = 10000  # 总用户数

# 创建用户特征模拟真实场景
df = pd.DataFrame({
    'user_id': range(n_users),
    'age': np.random.normal(35, 10, n_users).astype(int),
    'gender': np.random.choice(['M', 'F'], n_users, p=[0.55, 0.45]),
    'country': np.random.choice(['US', 'UK', 'CA', 'AU'], n_users, p=[0.6, 0.2, 0.1, 0.1]),
    'new_user': np.random.choice([0, 1], n_users, p=[0.3, 0.7])
})

# 分配用户到对照组A1和A2
df['group'] = np.random.choice(['A1', 'A2'], n_users, p=[0.5, 0.5])

# 模拟用户行为数据 - 两组应该相同
base_conversion = 0.15  # 基准转化率
# 添加一些噪声模拟真实环境
noise = np.random.normal(0, 0.01, n_users)

df['converted'] = 0
df.loc[df['new_user'] == 1, 'converted'] = np.random.binomial(1, base_conversion + 0.02 + noise[df['new_user'] == 1])
df.loc[df['new_user'] == 0, 'converted'] = np.random.binomial(1, base_conversion - 0.01 + noise[df['new_user'] == 0])

# 添加访问次数和停留时间
df['visit_count'] = np.random.poisson(3, n_users)
df['session_duration'] = np.random.gamma(5, 10, n_users)

print("数据前几行:")
print(df.head())
print(f"\n总用户数: {len(df)}")
print(f"A1组用户数: {len(df[df['group'] == 'A1'])}")
print(f"A2组用户数: {len(df[df['group'] == 'A2'])}")

# 检验分组平衡性
print("分组平衡性检验:")
print("=" * 50)

# 数值特征平衡检验
numeric_cols = ['age', 'visit_count', 'session_duration']
balance_results = []

for col in numeric_cols:
    group_a1 = df[df['group'] == 'A1'][col]
    group_a2 = df[df['group'] == 'A2'][col]
    
    t_stat, p_value = stats.ttest_ind(group_a1, group_a2)
    balance_results.append({
        'feature': col,
        'a1_mean': group_a1.mean(),
        'a2_mean': group_a2.mean(),
        'difference': group_a1.mean() - group_a2.mean(),
        'p_value': p_value
    })

balance_df = pd.DataFrame(balance_results)
print("数值特征平衡性:")
print(balance_df.round(4))

# 分类特征平衡检验
categorical_cols = ['gender', 'country', 'new_user']
cat_balance_results = []

for col in categorical_cols:
    cross_tab = pd.crosstab(df['group'], df[col])
    chi2, p_value, _, _ = stats.chi2_contingency(cross_tab)
    
    cat_balance_results.append({
        'feature': col,
        'chi2': chi2,
        'p_value': p_value
    })

cat_balance_df = pd.DataFrame(cat_balance_results)
print("\n分类特征平衡性:")
print(cat_balance_df.round(4))

# 可视化平衡性检查
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 数值特征分布
for i, col in enumerate(numeric_cols):
    sns.boxplot(x='group', y=col, data=df, ax=axes[0, i])
    axes[0, i].set_title(f'{col} 分布')

# 分类特征分布
for i, col in enumerate(categorical_cols[:3]):  # 只显示前3个分类变量
    prop_df = df.groupby(['group', col]).size().reset_index()
    prop_df['proportion'] = df.groupby(['group', col]).size().groupby(level=0).apply(lambda x: x / x.sum()).values
    sns.barplot(x='group', y='proportion', hue=col, data=prop_df, ax=axes[1, i])
    axes[1, i].set_title(f'{col} 比例分布')
    axes[1, i].legend_.remove()  # 移除图例避免拥挤

axes[1, 2].legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.savefig('balance_check.png', dpi=300, bbox_inches='tight')
plt.show()

# 核心指标分析 - 转化率
conversion_rates = df.groupby('group')['converted'].agg(['mean', 'count', 'std'])

# 计算置信区间
def get_confidence_interval(mean, std, count, confidence=0.95):
    se = std / np.sqrt(count)
    margin = se * stats.t.ppf((1 + confidence) / 2, count - 1)
    return mean - margin, mean + margin

conversion_rates['ci_lower'] = get_confidence_interval(
    conversion_rates['mean'], conversion_rates['std'], conversion_rates['count']
)[0]
conversion_rates['ci_upper'] = get_confidence_interval(
    conversion_rates['mean'], conversion_rates['std'], conversion_rates['count']
)[1]

print("转化率分析:")
print(conversion_rates.round(4))

# 统计显著性检验
a1_conversions = df[df['group'] == 'A1']['converted'].sum()
a1_total = len(df[df['group'] == 'A1'])
a2_conversions = df[df['group'] == 'A2']['converted'].sum()
a2_total = len(df[df['group'] == 'A2'])

z_stat, p_value = proportion.proportions_ztest(
    [a1_conversions, a2_conversions], 
    [a1_total, a2_total]
)

print(f"\n比例检验结果: z = {z_stat:.4f}, p = {p_value:.4f}")

# 可视化转化率比较
plt.figure(figsize=(10, 6))
sns.barplot(x='group', y='converted', data=df, ci=95)
plt.title('转化率比较 (A1 vs A2)')
plt.ylabel('转化率')
plt.xlabel('实验组')
plt.savefig('conversion_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

# 综合评估函数
def evaluate_aa_test(results, alpha=0.05):
    """
    综合评估A/A测试结果
    """
    # 检查所有平衡性检验
    balance_issues = len([r for r in results if r['p_value'] < alpha])
    
    # 检查主要指标显著性
    main_metric_sig = p_value < alpha
    
    # 评估系统状态
    if balance_issues > 0:
        return "FAIL", f"发现 {balance_issues} 个平衡性问题"
    elif main_metric_sig:
        return "FAIL", "主要指标显示显著差异"
    else:
        return "PASS", "系统表现正常"

# 收集所有检验结果
all_results = []
all_results.extend(balance_results)
all_results.extend(cat_balance_results.to_dict('records'))

# 执行评估
status, message = evaluate_aa_test(all_results)
print(f"\nA/A测试评估结果: {status}")
print(f"详细信息: {message}")

# 生成详细报告
print("\n" + "="*50)
print("A/A测试详细报告")
print("="*50)
print(f"样本总量: {n_users}")
print(f"分组比例: A1: {len(df[df['group'] == 'A1'])/n_users:.2%}, A2: {len(df[df['group'] == 'A2'])/n_users:.2%}")
print(f"总体转化率: {df['converted'].mean():.4f}")
print(f"转化率差异: {conversion_rates.loc['A1', 'mean'] - conversion_rates.loc['A2', 'mean']:.4f}")
print(f"统计显著性: {'是' if p_value < 0.05 else '否'} (p = {p_value:.4f})")
print(f"系统状态: {status} - {message}")

# 建议行动
if status == "PASS":
    print("\n建议: 实验系统可靠，可以开始A/B测试")
else:
    print("\n建议: 实验系统存在问题，需要先调查和修复")

# 样本量计算函数
def calculate_aa_sample_size(alpha=0.05, power=0.8, baseline_rate=0.15, mde=0.01):
    """
    计算A/A测试所需样本量
    alpha: 显著性水平
    power: 统计功效
    baseline_rate: 基准转化率
    mde: 最小可检测效应（绝对差异）
    """
    # 使用比例检验的样本量公式
    effect_size = proportion.proportion_effectsize(baseline_rate, baseline_rate + mde)
    sample_size = power.zt_ind_solve_power(
        effect_size=effect_size, 
        alpha=alpha, 
        power=power, 
        ratio=1.0
    )
    return int(sample_size)

# 示例计算
required_n = calculate_aa_sample_size(alpha=0.05, power=0.8, 
                                     baseline_rate=0.15, mde=0.01)
print(f"所需样本量(每组): {required_n}")
print(f"总样本量: {required_n * 2}")

# 不同参数下的样本量分析
parameters = [
    (0.05, 0.8, 0.10, 0.01),
    (0.05, 0.8, 0.15, 0.01),
    (0.05, 0.9, 0.15, 0.01),
    (0.01, 0.8, 0.15, 0.01)
]

print("\n不同参数下的样本量需求:")
print("α\t功效\t基准率\tMDE\t样本量")
for alpha, power, baseline, mde in parameters:
    n = calculate_aa_sample_size(alpha, power, baseline, mde)
    print(f"{alpha}\t{power}\t{baseline}\t{mde}\t{n}")

# 模拟多重检验问题
np.random.seed(42)
n_metrics = 20  # 检查20个指标
n_simulations = 1000  # 模拟1000次实验

# 模拟零假设为真时的p值分布
false_positives = []
for i in range(n_simulations):
    # 生成20个独立检验的p值
    p_values = [stats.ttest_ind(np.random.normal(0, 1, 100), 
                               np.random.normal(0, 1, 100)).pvalue 
               for _ in range(n_metrics)]
    # 记录是否有任何假阳性
    false_positives.append(any(p < 0.05 for p in p_values))

false_positive_rate = sum(false_positives) / n_simulations
print(f"未校正的假阳性率: {false_positive_rate:.3f}")

# 应用多重检验校正
def apply_multiple_testing_correction(p_values, method='bonferroni'):
    """
    应用多重检验校正
    """
    from statsmodels.stats.multitest import multipletests
    return multipletests(p_values, alpha=0.05, method=method)

# 比较不同校正方法
methods = ['bonferroni', 'holm', 'fdr_bh']
correction_results = {}

for method in methods:
    corrected_fp = 0
    for i in range(n_simulations):
        p_values = [stats.ttest_ind(np.random.normal(0, 1, 100), 
                                   np.random.normal(0, 1, 100)).pvalue 
                   for _ in range(n_metrics)]
        _, corrected, _, _ = apply_multiple_testing_correction(p_values, method)
        corrected_fp += any(corrected)
    
    correction_results[method] = corrected_fp / n_simulations

print("\n不同校正方法的假阳性率:")
for method, rate in correction_results.items():
    print(f"{method}: {rate:.3f}")

# 模拟他们发现的问题
np.random.seed(123)
n_users = 50000

# 模拟有bug的分流算法（基于用户ID的哈希函数问题）
def buggy_assignment(user_id):
    hash_val = hash(str(user_id)) % 100
    if hash_val < 30:  # bug: 应该是50，但实际只有30
        return 'A1'
    else:
        return 'A2'

user_ids = range(n_users)
groups = [buggy_assignment(uid) for uid in user_ids]

group_counts = pd.Series(groups).value_counts()
print("分流比例:")
print(group_counts / n_users)

# 模拟他们的发现
np.random.seed(42)
n_experiments = 100
p_values = []

for i in range(n_experiments):
    # 模拟两组完全相同的转化数据
    group_a = np.random.binomial(1, 0.12, 1000)
    group_b = np.random.binomial(1, 0.12, 1000)
    _, p_val = stats.ttest_ind(group_a, group_b)
    p_values.append(p_val)

# 分析p值分布
print(f"p < 0.05 的比例: {sum(np.array(p_values) < 0.05) / n_experiments:.3f}")
print(f"p值分布均匀性检验: {stats.kstest(p_values, 'uniform').pvalue:.4f}")

# 可视化
plt.figure(figsize=(10, 6))
plt.hist(p_values, bins=20, alpha=0.7, edgecolor='black')
plt.axhline(y=n_experiments/20, color='red', linestyle='--', label='期望值')
plt.xlabel('p值')
plt.ylabel频数')
plt.title('A/A测试中的p值分布')
plt.legend()
plt.show()