import numpy as np
from scipy.stats import norm

# 例: 工厂历史灯泡寿命 μ=1000h, σ=50h (已知)
# 抽 30 个, X̄ = 985h
n, xbar, mu0, sigma = 30, 985, 1000, 50
z = (xbar - mu0) / (sigma / np.sqrt(n))
p_two = 2 * (1 - norm.cdf(abs(z)))
print(f"z = {z:.3f}, p = {p_two:.4f}")
# z = -1.643, p = 0.100 → 不显著

from scipy.stats import t

# 同上例子, 但 σ 未知, 用 s 估计
n, xbar, mu0, s = 30, 985, 1000, 60  # s = 60 (从样本估计)
t_stat = (xbar - mu0) / (s / np.sqrt(n))
p_two = 2 * (1 - t.cdf(abs(t_stat), df=n-1))
print(f"t = {t_stat:.3f}, df = {n-1}, p = {p_two:.4f}")
# t = -1.369, p = 0.181 → 不显著 (因为 s=60 > σ=50, 噪声更大)

from scipy.stats import norm
n, p_hat, p0 = 500, 0.10, 0.08
z = (p_hat - p0) / np.sqrt(p0 * (1 - p0) / n)
p_two = 2 * (1 - norm.cdf(abs(z)))
print(f"z = {z:.3f}, p = {p_two:.4f}")
# z = 1.838, p = 0.066 → 临界, 0.05 水平下不显著

import numpy as np
from scipy import stats

# 1. 准备数据
np.random.seed(42)
sample = np.random.normal(98, 1.5, 25)  # 25 个测量值
mu0 = 100  # 标准值

# 2. 数据探索
print(f"n = {len(sample)}, X̄ = {sample.mean():.3f}, s = {sample.std(ddof=1):.3f}")
print(f"Min = {sample.min():.3f}, Max = {sample.max():.3f}")

# 3. 检验正态性 (小样本必做)
stat, p_norm = stats.shapiro(sample)
print(f"Shapiro-Wilk 正态性: stat = {stat:.3f}, p = {p_norm:.3f}")
# p > 0.05 → 不拒绝正态性, 可用 t 检验

# 4. t 检验 (scipy 一行)
t_stat, p_two = stats.ttest_1samp(sample, mu0)
print(f"t 检验: t = {t_stat:.3f}, p = {p_two:.4f}")

# 5. 决策
alpha = 0.05
if p_two < alpha:
    print(f"✅ p < {alpha}, 拒绝 H₀: 样本均值显著不等于 {mu0}")
else:
    print(f"❌ p ≥ {alpha}, 不能拒绝 H₀")

# 6. 95% 置信区间 (配套)
t_crit = stats.t.ppf(0.975, len(sample) - 1)
margin = t_crit * sample.std(ddof=1) / np.sqrt(len(sample))
print(f"95% CI: [{sample.mean() - margin:.2f}, {sample.mean() + margin:.2f}]")

from statsmodels.stats.power import ttest_power
# 已知: 效应大小 Cohen's d = 1/2 = 0.5
# 找 n
from statsmodels.stats.power import TTestIndPower
analysis = TTestIndPower()
n_required = analysis.solve_power(effect_size=0.5, alpha=0.05, power=0.8)
print(f"达到 80% power 所需 n = {n_required:.0f}")  # ≈ 64

# 反例: 显著但效应极小
n = 100000
sample_a = np.random.normal(0, 1, n)
sample_b = np.random.normal(0.01, 1, n)  # 差 0.01
t_stat, p_val = stats.ttest_ind(sample_a, sample_b)
print(f"t = {t_stat:.3f}, p = {p_val:.4f}")
# p ≈ 0.04 (显著), 但效应大小 0.01/1 = 0.01, 完全不重要