概率与统计（七）：假设检验——p 值、置信区间及其全部陷阱

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(42)

data = np.random.normal(loc=10.2, scale=0.5, size=25)
t_stat, p_value = stats.ttest_1samp(data, popmean=10.0)
print(f"One-sample t-test: t = {t_stat:.4f}, p = {p_value:.4f}")

group_a = np.random.normal(loc=5.0, scale=1.0, size=30)
group_b = np.random.normal(loc=5.5, scale=1.2, size=35)
t_stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"Two-sample t-test: t = {t_stat:.4f}, p = {p_value:.4f}")

observed = np.array([[30, 10], [15, 25]])
chi2, p_value, dof, expected = stats.chi2_contingency(observed)
print(f"Chi-squared test: chi2 = {chi2:.4f}, p = {p_value:.4f}, dof = {dof}")

n_A, n_B = 5000, 5000
conversions_A, conversions_B = 500, 560
p_A = conversions_A / n_A
p_B = conversions_B / n_B
p_pool = (conversions_A + conversions_B) / (n_A + n_B)
se = np.sqrt(p_pool * (1 - p_pool) * (1/n_A + 1/n_B))
z = (p_B - p_A) / se
p_value = 2 * (1 - stats.norm.cdf(abs(z)))
print(f"A/B test: z = {z:.4f}, p = {p_value:.4f}")
print(f"  Lift: {p_B - p_A:.4f} ({(p_B - p_A)/p_A*100:.1f}%)")
print(f"  95% CI for lift: [{p_B-p_A-1.96*se:.4f}, {p_B-p_A+1.96*se:.4f}]")

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

ax = axes[0]
x = np.linspace(-4, 4, 300)
ax.plot(x, stats.norm.pdf(x), 'b-', linewidth=2)
t_obs = 2.0
ax.fill_between(x[x >= t_obs], stats.norm.pdf(x[x >= t_obs]),
                alpha=0.4, color='red', label=f'p-value/2 (right tail)')
ax.fill_between(x[x <= -t_obs], stats.norm.pdf(x[x <= -t_obs]),
                alpha=0.4, color='red', label=f'p-value/2 (left tail)')
ax.axvline(t_obs, color='red', linestyle='--', alpha=0.7)
ax.axvline(-t_obs, color='red', linestyle='--', alpha=0.7)
ax.set_title(f'Two-sided p-value (z = {t_obs})', fontsize=13)
ax.set_xlabel('Test statistic')
ax.set_ylabel('Density under $H_0$')
ax.legend(fontsize=10)

ax = axes[1]
x = np.linspace(-4, 6, 300)
ax.plot(x, stats.norm.pdf(x, 0, 1), 'b-', linewidth=2, label='$H_0$: $\\mu = 0$')
delta = 2.0
ax.plot(x, stats.norm.pdf(x, delta, 1), 'r-', linewidth=2,
        label=f'$H_1$: $\\mu = {delta}$')
z_crit = 1.645  # one-sided alpha = 0.05
ax.axvline(z_crit, color='gray', linestyle='--')
ax.fill_between(x[x >= z_crit], stats.norm.pdf(x[x >= z_crit], delta, 1),
                alpha=0.3, color='green', label=f'Power = {1-stats.norm.cdf(z_crit-delta):.3f}')
ax.fill_between(x[x >= z_crit], stats.norm.pdf(x[x >= z_crit], 0, 1),
                alpha=0.3, color='red', label=f'$\\alpha$ = {1-stats.norm.cdf(z_crit):.3f}')
ax.set_title('Power of a Test', fontsize=13)
ax.set_xlabel('Test statistic')
ax.legend(fontsize=9)

ax = axes[2]
m_tests = np.arange(1, 101)
fwer = 1 - (1 - 0.05)**m_tests
ax.plot(m_tests, fwer, 'r-', linewidth=2, label='FWER (no correction)')
ax.axhline(0.05, color='gray', linestyle=':', label='$\\alpha$ = 0.05')
ax.set_xlabel('Number of tests', fontsize=12)
ax.set_ylabel('P(at least one false positive)', fontsize=12)
ax.set_title('Multiple Testing Problem', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hypothesis_testing.png', dpi=150)
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

from scipy.stats import false_discovery_control

np.random.seed(42)
n_tests = 100
n_true_effects = 10
n = 30  # samples per test

p_values = []
truth = []

for i in range(n_tests):
    if i < n_true_effects:
        data = np.random.normal(loc=2.0, scale=1, size=n)
        truth.append(True)
    else:
        data = np.random.normal(loc=0.0, scale=1, size=n)
        truth.append(False)
    _, pval = stats.ttest_1samp(data, 0)
    p_values.append(pval)

p_values = np.array(p_values)
truth = np.array(truth)

bonferroni_reject = p_values < 0.05 / n_tests

sorted_idx = np.argsort(p_values)
bh_threshold = np.arange(1, n_tests + 1) / n_tests * 0.05
bh_reject_sorted = p_values[sorted_idx] <= bh_threshold
if bh_reject_sorted.any():
    k = np.max(np.where(bh_reject_sorted)[0]) + 1
    bh_reject = np.zeros(n_tests, dtype=bool)
    bh_reject[sorted_idx[:k]] = True
else:
    bh_reject = np.zeros(n_tests, dtype=bool)

print(f"{'Method':<25} {'Rejections':>10} {'True Pos':>10} {'False Pos':>10}")
print("-" * 55)
print(f"{'No correction':<25} {(p_values < 0.05).sum():>10} "
      f"{((p_values < 0.05) & truth).sum():>10} "
      f"{((p_values < 0.05) & ~truth).sum():>10}")
print(f"{'Bonferroni':<25} {bonferroni_reject.sum():>10} "
      f"{(bonferroni_reject & truth).sum():>10} "
      f"{(bonferroni_reject & ~truth).sum():>10}")
print(f"{'Benjamini-Hochberg':<25} {bh_reject.sum():>10} "
      f"{(bh_reject & truth).sum():>10} "
      f"{(bh_reject & ~truth).sum():>10}")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import numpy as np
from scipy import stats

np.random.seed(42)

group_a = np.array([5.1, 4.8, 5.3, 4.9, 5.0, 5.2, 5.4, 4.7])
group_b = np.array([5.5, 5.8, 5.3, 5.7, 5.9, 5.6, 5.4, 5.8])

observed_diff = group_b.mean() - group_a.mean()
combined = np.concatenate([group_a, group_b])
n_a = len(group_a)
n_perms = 100_000

perm_diffs = np.zeros(n_perms)
for i in range(n_perms):
    np.random.shuffle(combined)
    perm_diffs[i] = combined[n_a:].mean() - combined[:n_a].mean()

p_value = np.mean(np.abs(perm_diffs) >= np.abs(observed_diff))
print(f"Observed difference: {observed_diff:.4f}")
print(f"Permutation p-value: {p_value:.4f}")

t_stat, t_pval = stats.ttest_ind(group_a, group_b)
print(f"t-test p-value: {t_pval:.4f}")