统计学基础与假设检验

import numpy as np
import pandas as pd

data = [23, 25, 27, 23, 29, 31, 23, 35, 28, 26]

# 均值（Mean）- 算术平均
mean = np.mean(data)
print(f"均值: {mean:.2f}")  # 27.00

# 中位数（Median）- 排序后中间的值
median = np.median(data)
print(f"中位数: {median:.2f}")  # 26.50

# 众数（Mode）- 出现最频繁的值
from scipy import stats
mode = stats.mode(data, keepdims=True)
print(f"众数: {mode.mode[0]}")  # 23

# 方差（Variance）- 数据与均值的平方差的平均
variance = np.var(data, ddof=1)  # ddof=1 表示样本方差
print(f"方差: {variance:.2f}")

# 标准差（Standard Deviation）- 方差的平方根
std = np.std(data, ddof=1)
print(f"标准差: {std:.2f}")

# 四分位距（IQR）- 抗极端值
Q1 = np.percentile(data, 25)
Q3 = np.percentile(data, 75)
IQR = Q3 - Q1
print(f"四分位距: {IQR:.2f}")

# 变异系数（Coefficient of Variation）- 标准化的离散度
CV = (std / mean) * 100
print(f"变异系数: {CV:.2f}%")

# 偏度（Skewness）- 衡量对称性
skewness = stats.skew(data)
print(f"偏度: {skewness:.2f}")
# < 0: 左偏（负偏），> 0: 右偏（正偏）

# 峰度（Kurtosis）- 衡量尾部厚度
kurtosis = stats.kurtosis(data)
print(f"峰度: {kurtosis:.2f}")
# > 0: 尖峰，< 0: 平峰

import matplotlib.pyplot as plt
import seaborn as sns

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 直方图 + KDE
sns.histplot(data, kde=True, ax=axes[0,0])
axes[0,0].set_title('分布直方图')

# 箱线图（Box Plot）
sns.boxplot(data=data, ax=axes[0,1])
axes[0,1].set_title('箱线图')

# Q-Q图（检验正态性）
stats.probplot(data, dist="norm", plot=axes[1,0])
axes[1,0].set_title('Q-Q图')

# 小提琴图（Violin Plot）
sns.violinplot(data=data, ax=axes[1,1])
axes[1,1].set_title('小提琴图')

plt.tight_layout()
plt.show()

# 模拟中心极限定理
population = np.random.exponential(scale=2, size=100000)  # 总体：指数分布

# 重复抽样
sample_means = []
for _ in range(1000):
    sample = np.random.choice(population, size=30)
    sample_means.append(np.mean(sample))

# 样本均值的分布趋向正态分布
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.hist(population, bins=50, alpha=0.7, edgecolor='black')
plt.title('总体分布（指数分布）')
plt.xlabel('值')
plt.ylabel('频数')

plt.subplot(1, 2, 2)
plt.hist(sample_means, bins=30, alpha=0.7, edgecolor='black')
plt.title('样本均值分布（趋向正态）')
plt.xlabel('样本均值')
plt.ylabel('频数')

plt.tight_layout()
plt.show()

print(f"总体均值: {np.mean(population):.2f}")
print(f"样本均值的均值: {np.mean(sample_means):.2f}")
print(f"样本均值的标准误: {np.std(sample_means):.2f}")

from scipy import stats

def confidence_interval(data, confidence=0.95):
    """计算均值的置信区间"""
    n = len(data)
    mean = np.mean(data)
    std_err = stats.sem(data)  # 标准误 = std / sqrt(n)
    
    # t分布（样本量小时）或正态分布
    if n < 30:
        margin_error = std_err * stats.t.ppf((1 + confidence) / 2, n - 1)
    else:
        margin_error = std_err * stats.norm.ppf((1 + confidence) / 2)
    
    ci_lower = mean - margin_error
    ci_upper = mean + margin_error
    
    return mean, (ci_lower, ci_upper)

# 示例
data = [23, 25, 27, 23, 29, 31, 23, 35, 28, 26]
mean, (ci_lower, ci_upper) = confidence_interval(data, confidence=0.95)

print(f"样本均值: {mean:.2f}")
print(f"95% 置信区间: [{ci_lower:.2f}, {ci_upper:.2f}]")
# 解释：我们有95%的把握认为总体均值在此区间内

from scipy import stats

# 实测数据
volumes = [498, 502, 496, 501, 499, 497, 503, 500, 498, 502]

# H₀: μ = 500（实际容量符合标称）
# H₁: μ < 500（实际容量少于标称）

claimed_mean = 500

# 单侧t检验
t_statistic, p_value = stats.ttest_1samp(volumes, claimed_mean)

# 单侧p值（仅关心"少于"的情况）
p_value_one_sided = p_value / 2 if t_statistic < 0 else 1 - p_value / 2

print(f"样本均值: {np.mean(volumes):.2f} ml")
print(f"t统计量: {t_statistic:.4f}")
print(f"p值（单侧）: {p_value_one_sided:.4f}")

if p_value_one_sided < 0.05:
    print("拒绝原假设：存在缺斤短两问题")
else:
    print("不能拒绝原假设：未发现缺斤短两")

# Cohen's d：标准化的均值差异
cohens_d = (np.mean(volumes) - claimed_mean) / np.std(volumes, ddof=1)
print(f"Cohen's d: {cohens_d:.4f}")

# 解释：
# |d| < 0.2: 小效应
# 0.2 ≤ |d| < 0.8: 中效应
# |d| ≥ 0.8: 大效应

# 策略A：传统广告
strategy_a = [120, 135, 128, 142, 130, 125, 138, 132, 140, 127]

# 策略B：创新广告
strategy_b = [145, 152, 138, 160, 148, 155, 142, 158, 150, 147]

# H₀: μ_A = μ_B（两种策略效果相同）
# H₁: μ_A ≠ μ_B（两种策略效果不同）

# 先检验方差齐性（Levene's Test）
stat, p_levene = stats.levene(strategy_a, strategy_b)
print(f"Levene检验p值: {p_levene:.4f}")

if p_levene > 0.05:
    # 方差齐性，使用标准t检验
    t_stat, p_value = stats.ttest_ind(strategy_a, strategy_b)
    print("使用标准t检验（方差齐性）")
else:
    # 方差不齐，使用Welch's t检验
    t_stat, p_value = stats.ttest_ind(strategy_a, strategy_b, equal_var=False)
    print("使用Welch's t检验（方差不齐）")

print(f"t统计量: {t_stat:.4f}")
print(f"p值: {p_value:.4f}")

# 效应量
pooled_std = np.sqrt((np.var(strategy_a, ddof=1) + np.var(strategy_b, ddof=1)) / 2)
cohens_d = (np.mean(strategy_b) - np.mean(strategy_a)) / pooled_std
print(f"Cohen's d: {cohens_d:.4f}")

# 可视化
plt.figure(figsize=(10, 6))
plt.boxplot([strategy_a, strategy_b], labels=['策略A', '策略B'])
plt.ylabel('转化量')
plt.title('两种广告策略效果对比')
plt.grid(axis='y', alpha=0.3)
plt.show()

if p_value < 0.05:
    print(f"结论：策略B显著优于策略A（p={p_value:.4f}）")
else:
    print("结论：两种策略无显著差异")

# 培训前绩效
before = [72, 68, 75, 70, 73, 69, 71, 74, 68, 72]

# 培训后绩效（同一批员工）
after = [78, 72, 80, 75, 79, 73, 76, 80, 74, 77]

# H₀: μ_after - μ_before = 0（培训无效）
# H₁: μ_after - μ_before > 0（培训有效）

t_stat, p_value = stats.ttest_rel(after, before)

# 单侧p值
p_value_one_sided = p_value / 2

print(f"平均提升: {np.mean(after) - np.mean(before):.2f} 分")
print(f"t统计量: {t_stat:.4f}")
print(f"p值（单侧）: {p_value_one_sided:.4f}")

# 可视化差异
differences = np.array(after) - np.array(before)

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
x = range(len(before))
plt.plot(x, before, 'o-', label='培训前', markersize=8)
plt.plot(x, after, 's-', label='培训后', markersize=8)
plt.xlabel('员工编号')
plt.ylabel('绩效得分')
plt.legend()
plt.title('培训前后绩效对比')
plt.grid(alpha=0.3)

plt.subplot(1, 2, 2)
plt.bar(range(len(differences)), differences, color='steelblue', alpha=0.7)
plt.axhline(0, color='red', linestyle='--', linewidth=2)
plt.xlabel('员工编号')
plt.ylabel('绩效提升')
plt.title('每个员工的绩效变化')
plt.grid(axis='y', alpha=0.3)

plt.tight_layout()
plt.show()

if p_value_one_sided < 0.05:
    print("结论：培训显著提升了绩效")
else:
    print("结论：培训效果不显著")

# 三种教学方法的测试成绩
method_a = [85, 88, 82, 90, 87, 84, 89, 86]
method_b = [78, 82, 75, 80, 79, 77, 81, 76]
method_c = [92, 95, 90, 93, 94, 91, 96, 93]

# H₀: μ_A = μ_B = μ_C（三种方法效果相同）
# H₁: 至少有两组均值不同

# 单因素方差分析
f_stat, p_value = stats.f_oneway(method_a, method_b, method_c)

print(f"F统计量: {f_stat:.4f}")
print(f"p值: {p_value:.4f}")

if p_value < 0.05:
    print("结论：至少有两种方法的效果存在显著差异")
    
    # 事后多重比较（Tukey HSD）
    from scipy.stats import tukey_hsd
    
    all_data = np.concatenate([method_a, method_b, method_c])
    groups = ['A']*len(method_a) + ['B']*len(method_b) + ['C']*len(method_c)
    
    # 使用statsmodels进行多重比较
    from statsmodels.stats.multicomp import pairwise_tukeyhsd
    
    tukey_result = pairwise_tukeyhsd(all_data, groups, alpha=0.05)
    print("\nTukey HSD事后比较：")
    print(tukey_result)
else:
    print("结论：三种方法效果无显著差异")

# 可视化
plt.figure(figsize=(10, 6))
plt.boxplot([method_a, method_b, method_c], labels=['方法A', '方法B', '方法C'])
plt.ylabel('测试成绩')
plt.title('三种教学方法效果对比')
plt.grid(axis='y', alpha=0.3)
plt.show()

def eta_squared(groups):
    """计算效应量η²"""
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)
    
    # 组间平方和（SS_between）
    ss_between = sum(len(group) * (np.mean(group) - grand_mean)**2 
                     for group in groups)
    
    # 总平方和（SS_total）
    ss_total = sum((x - grand_mean)**2 for x in all_data)
    
    # η² = SS_between / SS_total
    eta_sq = ss_between / ss_total
    
    return eta_sq

eta_sq = eta_squared([method_a, method_b, method_c])
print(f"效应量 η²: {eta_sq:.4f}")

# 解释：
# η² < 0.01: 小效应
# 0.01 ≤ η² < 0.06: 中效应
# η² ≥ 0.06: 大效应

# 投掷600次的结果
observed = [95, 105, 98, 102, 110, 90]  # 各点数出现次数

# H₀: 骰子是公平的（各点数概率均为1/6）
expected = [100, 100, 100, 100, 100, 100]  # 理论频数

chi2_stat, p_value = stats.chisquare(observed, expected)

print(f"卡方统计量: {chi2_stat:.4f}")
print(f"p值: {p_value:.4f}")

if p_value < 0.05:
    print("结论：骰子不公平")
else:
    print("结论：骰子是公平的")

# 可视化
x = range(1, 7)
width = 0.35

plt.figure(figsize=(10, 6))
plt.bar([i - width/2 for i in x], observed, width, label='观测频数', alpha=0.8)
plt.bar([i + width/2 for i in x], expected, width, label='期望频数', alpha=0.8)
plt.xlabel('骰子点数')
plt.ylabel('出现次数')
plt.title('骰子投掷结果分析')
plt.legend()
plt.xticks(x)
plt.grid(axis='y', alpha=0.3)
plt.show()

# 列联表（Contingency Table）
#           产品A  产品B  产品C
# 男性        30    20    50
# 女性        45    35    20

observed = np.array([
    [30, 20, 50],  # 男性
    [45, 35, 20]   # 女性
])

# H₀: 性别与产品偏好独立
# H₁: 性别与产品偏好相关

chi2_stat, p_value, dof, expected = stats.chi2_contingency(observed)

print(f"卡方统计量: {chi2_stat:.4f}")
print(f"自由度: {dof}")
print(f"p值: {p_value:.4f}")
print(f"\n期望频数：\n{expected}")

if p_value < 0.05:
    print("\n结论：性别与产品偏好相关")
else:
    print("\n结论：性别与产品偏好独立")

# 效应量：Cramér's V
n = observed.sum()
min_dim = min(observed.shape) - 1
cramers_v = np.sqrt(chi2_stat / (n * min_dim))
print(f"Cramér's V: {cramers_v:.4f}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 观测频数
pd.DataFrame(observed, 
             index=['男性', '女性'],
             columns=['产品A', '产品B', '产品C']).plot(kind='bar', ax=axes[0], alpha=0.8)
axes[0].set_title('观测频数')
axes[0].set_ylabel('人数')
axes[0].legend(title='产品')
axes[0].grid(axis='y', alpha=0.3)

# 期望频数
pd.DataFrame(expected, 
             index=['男性', '女性'],
             columns=['产品A', '产品B', '产品C']).plot(kind='bar', ax=axes[1], alpha=0.8)
axes[1].set_title('期望频数（独立假设下）')
axes[1].set_ylabel('人数')
axes[1].legend(title='产品')
axes[1].grid(axis='y', alpha=0.3)

plt.tight_layout()
plt.show()

from statsmodels.stats.power import tt_ind_solve_power

def calculate_sample_size(effect_size, alpha=0.05, power=0.8):
    """
    计算所需样本量
    
    参数：
    - effect_size: Cohen's d（效应量）
    - alpha: 显著性水平（第一类错误率）
    - power: 检验功效（1 - 第二类错误率）
    """
    n = tt_ind_solve_power(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        alternative='two-sided'
    )
    
    return int(np.ceil(n))

# 示例：检测5%的提升（假设基准转化率10%）
baseline_rate = 0.10
target_rate = 0.105  # 提升5%

# 计算效应量
pooled_std = np.sqrt((baseline_rate * (1 - baseline_rate) + 
                       target_rate * (1 - target_rate)) / 2)
effect_size = (target_rate - baseline_rate) / pooled_std

n_per_group = calculate_sample_size(effect_size, alpha=0.05, power=0.8)
print(f"效应量: {effect_size:.4f}")
print(f"每组所需样本量: {n_per_group}")
print(f"总样本量: {n_per_group * 2}")

class ABTest:
    def __init__(self, control_data, treatment_data, metric_type='continuous'):
        """
        A/B测试分析
        
        参数：
        - control_data: 对照组数据
        - treatment_data: 实验组数据
        - metric_type: 'continuous'（连续）或 'conversion'（转化率）
        """
        self.control = np.array(control_data)
        self.treatment = np.array(treatment_data)
        self.metric_type = metric_type
        
    def run_test(self, alpha=0.05):
        """执行假设检验"""
        if self.metric_type == 'continuous':
            return self._test_continuous(alpha)
        elif self.metric_type == 'conversion':
            return self._test_conversion(alpha)
    
    def _test_continuous(self, alpha):
        """连续指标检验（如收入、时长）"""
        # t检验
        t_stat, p_value = stats.ttest_ind(self.treatment, self.control)
        
        # 效应量
        pooled_std = np.sqrt((np.var(self.control, ddof=1) + 
                              np.var(self.treatment, ddof=1)) / 2)
        cohens_d = (np.mean(self.treatment) - np.mean(self.control)) / pooled_std
        
        # 置信区间
        mean_diff = np.mean(self.treatment) - np.mean(self.control)
        se_diff = np.sqrt(np.var(self.control, ddof=1)/len(self.control) + 
                          np.var(self.treatment, ddof=1)/len(self.treatment))
        ci_lower = mean_diff - 1.96 * se_diff
        ci_upper = mean_diff + 1.96 * se_diff
        
        return {
            'test_type': 't-test',
            'control_mean': np.mean(self.control),
            'treatment_mean': np.mean(self.treatment),
            'mean_difference': mean_diff,
            'relative_lift': (mean_diff / np.mean(self.control)) * 100,
            'p_value': p_value,
            'significant': p_value < alpha,
            'cohens_d': cohens_d,
            'confidence_interval': (ci_lower, ci_upper)
        }
    
    def _test_conversion(self, alpha):
        """转化率指标检验"""
        # 转化数
        control_conv = np.sum(self.control)
        treatment_conv = np.sum(self.treatment)
        
        # 总数
        control_n = len(self.control)
        treatment_n = len(self.treatment)
        
        # 转化率
        control_rate = control_conv / control_n
        treatment_rate = treatment_conv / treatment_n
        
        # z检验（比例检验）
        pooled_rate = (control_conv + treatment_conv) / (control_n + treatment_n)
        se = np.sqrt(pooled_rate * (1 - pooled_rate) * (1/control_n + 1/treatment_n))
        z_stat = (treatment_rate - control_rate) / se
        p_value = 2 * (1 - stats.norm.cdf(abs(z_stat)))
        
        # 置信区间
        se_diff = np.sqrt(control_rate*(1-control_rate)/control_n + 
                          treatment_rate*(1-treatment_rate)/treatment_n)
        ci_lower = (treatment_rate - control_rate) - 1.96 * se_diff
        ci_upper = (treatment_rate - control_rate) + 1.96 * se_diff
        
        return {
            'test_type': 'z-test (proportion)',
            'control_rate': control_rate,
            'treatment_rate': treatment_rate,
            'absolute_lift': treatment_rate - control_rate,
            'relative_lift': ((treatment_rate - control_rate) / control_rate) * 100,
            'p_value': p_value,
            'significant': p_value < alpha,
            'confidence_interval': (ci_lower, ci_upper)
        }
    
    def visualize(self):
        """可视化结果"""
        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        
        # 分布对比
        if self.metric_type == 'continuous':
            axes[0].hist(self.control, alpha=0.6, label='对照组', bins=20)
            axes[0].hist(self.treatment, alpha=0.6, label='实验组', bins=20)
            axes[0].axvline(np.mean(self.control), color='blue', linestyle='--', 
                           label=f'对照组均值: {np.mean(self.control):.2f}')
            axes[0].axvline(np.mean(self.treatment), color='orange', linestyle='--', 
                           label=f'实验组均值: {np.mean(self.treatment):.2f}')
        else:
            rates = [np.mean(self.control), np.mean(self.treatment)]
            axes[0].bar(['对照组', '实验组'], rates, alpha=0.7)
            for i, rate in enumerate(rates):
                axes[0].text(i, rate, f'{rate:.2%}', ha='center', va='bottom')
        
        axes[0].set_title('数据分布对比')
        axes[0].legend()
        axes[0].grid(axis='y', alpha=0.3)
        
        # 箱线图
        axes[1].boxplot([self.control, self.treatment], labels=['对照组', '实验组'])
        axes[1].set_title('箱线图对比')
        axes[1].grid(axis='y', alpha=0.3)
        
        plt.tight_layout()
        plt.show()

# 使用示例：连续指标
control_revenue = np.random.normal(100, 20, 1000)  # 对照组收入
treatment_revenue = np.random.normal(105, 20, 1000)  # 实验组收入（提升5元）

ab_test = ABTest(control_revenue, treatment_revenue, metric_type='continuous')
results = ab_test.run_test(alpha=0.05)

print("A/B测试结果：")
print(f"对照组均值: {results['control_mean']:.2f}")
print(f"实验组均值: {results['treatment_mean']:.2f}")
print(f"绝对提升: {results['mean_difference']:.2f}")
print(f"相对提升: {results['relative_lift']:.2f}%")
print(f"p值: {results['p_value']:.4f}")
print(f"显著性: {'显著' if results['significant'] else '不显著'}")
print(f"效应量(Cohen's d): {results['cohens_d']:.4f}")
print(f"95%置信区间: [{results['confidence_interval'][0]:.2f}, {results['confidence_interval'][1]:.2f}]")

ab_test.visualize()

# 使用示例：转化率指标
control_conversions = np.random.binomial(1, 0.10, 10000)  # 10%转化率
treatment_conversions = np.random.binomial(1, 0.105, 10000)  # 10.5%转化率

ab_test_conv = ABTest(control_conversions, treatment_conversions, metric_type='conversion')
results_conv = ab_test_conv.run_test(alpha=0.05)

print("\n转化率A/B测试结果：")
print(f"对照组转化率: {results_conv['control_rate']:.2%}")
print(f"实验组转化率: {results_conv['treatment_rate']:.2%}")
print(f"绝对提升: {results_conv['absolute_lift']:.2%}")
print(f"相对提升: {results_conv['relative_lift']:.2f}%")
print(f"p值: {results_conv['p_value']:.4f}")
print(f"显著性: {'显著' if results_conv['significant'] else '不显著'}")

def bonferroni_correction(p_values, alpha=0.05):
    """Bonferroni校正"""
    n_tests = len(p_values)
    adjusted_alpha = alpha / n_tests
    
    results = []
    for i, p in enumerate(p_values):
        results.append({
            'test': i + 1,
            'p_value': p,
            'adjusted_alpha': adjusted_alpha,
            'significant': p < adjusted_alpha
        })
    
    return results

# 示例：3个A/B测试
p_values = [0.03, 0.02, 0.045]

print("Bonferroni校正前：")
for i, p in enumerate(p_values):
    print(f"测试{i+1}: p={p:.4f}, 显著: {p < 0.05}")

print("\nBonferroni校正后：")
corrected_results = bonferroni_correction(p_values, alpha=0.05)
for result in corrected_results:
    print(f"测试{result['test']}: p={result['p_value']:.4f}, "
          f"调整后α={result['adjusted_alpha']:.4f}, "
          f"显著: {result['significant']}")

# 示例：大样本量下的"显著但无意义"结果
control_large = np.random.normal(100, 10, 100000)
treatment_large = np.random.normal(100.5, 10, 100000)  # 仅提升0.5

t_stat, p_value = stats.ttest_ind(treatment_large, control_large)
cohens_d = (np.mean(treatment_large) - np.mean(control_large)) / np.std(control_large)

print(f"p值: {p_value:.6f} (显著)")
print(f"效应量: {cohens_d:.6f} (极小)")
print(f"结论: 统计显著但实际意义很小")

# 检验功效分析
from statsmodels.stats.power import tt_ind_solve_power

effect_sizes = [0.2, 0.5, 0.8]  # 小、中、大效应
sample_sizes = [10, 30, 50, 100, 200]

plt.figure(figsize=(10, 6))

for effect_size in effect_sizes:
    powers = []
    for n in sample_sizes:
        power = tt_ind_solve_power(
            effect_size=effect_size,
            nobs1=n,
            alpha=0.05,
            alternative='two-sided'
        )
        powers.append(power)
    
    plt.plot(sample_sizes, powers, marker='o', 
             label=f'效应量 d={effect_size}')

plt.axhline(0.8, color='red', linestyle='--', label='目标功效 (80%)')
plt.xlabel('样本量')
plt.ylabel('检验功效 (Power)')
plt.title('样本量与检验功效关系')
plt.legend()
plt.grid(alpha=0.3)
plt.show()