9.3 假設檢定的計算實作#

import numpy as np, math

rng = np.random.default_rng(0)

# --- t 檢定：完整手作實作 ---
def t_test(a, b):
    n1, n2 = len(a), len(b)
    m1, m2 = a.mean(), b.mean()
    v1, v2 = a.var(ddof=1), b.var(ddof=1)
    se = math.sqrt(v1/n1 + v2/n2)
    t = (m1 - m2) / se
    # Welch 自由度
    df = (v1/n1 + v2/n2)**2 / ((v1/n1)**2/(n1-1) + (v2/n2)**2/(n2-1))
    # t 分布尾機率：用不完全貝塔函數
    from math import lgamma
    def t_sf(x, df):
        # 右尾機率，透過貝塔函數恆等式
        xx = df / (df + x*x)
        ib = _ibeta(xx, df/2, 0.5)
        return ib
    def _ibeta(x, a, b):
        # 正則化不完全貝塔函數（Lentz 連分數，簡化版）
        if x > (a + 1) / (a + b + 2):
            return 1 - _ibeta(1 - x, b, a)
        lbeta = lgamma(a+b) - lgamma(a) - lgamma(b)
        front = math.exp(a*math.log(x) + lbeta) / a
        f, c, d = 1.0, 0.0, 0.0
        for i in range(200):
            m = i // 2
            num = m*(b-m)*x / ((a+2*m-1)*(a+2*m)) if i % 2 == 0 \
                  else -((a+m)*(a+b+m)*x) / ((a+2*m)*(a+2*m+1))
            d = 1 + num*d
            if abs(d) < 1e-30: d = 1e-30
            d = 1 / d
            c = 1 + num/c
            if abs(c) < 1e-30: c = 1e-30
            f *= c * d
            if abs(c*d - 1) < 1e-12: break
        return min(max(front * f, 0), 1)
    p = 2 * t_sf(abs(t), df)
    return t, df, p

a = rng.normal(loc=100, scale=15, size=40)     # 對照組
b = rng.normal(loc=106, scale=15, size=45)     # 實驗組
t, df, p = t_test(a, b)
print(f"t 檢定: t = {t:.3f}, df = {df:.1f}, p = {p:.4f}")
try:
    from scipy import stats
    ref = stats.ttest_ind(a, b, equal_var=False)
    print(f"scipy 對照: t = {ref.statistic:.3f}, p = {ref.pvalue:.4f}")
except ImportError:
    print("（未安裝 scipy，略過對照）")

# --- 卡方檢定：骰子是否公正 ---
obs = np.bincount(rng.integers(1, 7, size=600), minlength=7)[1:]
exp = np.full(6, 100.0)
chi2 = ((obs - exp)**2 / exp).sum()
def chi2_sf(x, k):
    """卡方右尾機率（正則化不完全伽瑪 Q(k/2, x/2)）"""
    from math import lgamma
    s, term = 0.0, 1.0
    a, xx = k/2, x/2
    # 用級數算 P(a, xx)，Q = 1 - P（x 適中時夠用）
    if xx < a + 1:
        term = 1.0 / a; s = term; n = a
        for _ in range(500):
            n += 1; term *= xx / n; s += term
            if abs(term) < 1e-15 * abs(s): break
        return 1 - math.exp(-xx + a*math.log(xx) - lgamma(a)) * s
    return 1.0  # 簡化：此例 x 不會太大
print(f"骰子觀察 {obs.tolist()}，χ² = {chi2:.2f}, p ≈ {chi2_sf(chi2, 5):.4f}")

# --- 置換檢定：無分布假設的 p 值 ---
def permutation_test(a, b, n_perm=100_000):
    obs_diff = abs(a.mean() - b.mean())
    pool = np.concatenate([a, b])
    count = 0
    for _ in range(n_perm):
        rng.shuffle(pool)
        if abs(pool[:len(a)].mean() - pool[len(a):].mean()) >= obs_diff:
            count += 1
    return count / n_perm

print(f"置換檢定 p 值 ≈ {permutation_test(a, b):.4f}（應與 t 檢定接近）")