일표본 평균 검정¶

import math
from scipy.stats import t as tdist, norm

def test_mean_one_sample(xbar, n, mu0=0.0, sd=None, known_sigma=None,
                         alternative="two-sided", alpha=0.05):
    """known_sigma를 주면 z-검정, 아니면 표본 sd로 t-검정.

    원자료가 아니라 요약통계량(xbar, n, sd)만 받는다.
    검정에 필요한 것이 그것뿐이기 때문이다.
    alternative는 scipy의 관례를 그대로 따른다.
    돌려주는 값은 (통계량, p-값, 기각 여부, 이름)이다.
    """
    if known_sigma is not None:
        se = known_sigma / math.sqrt(n)
        z = (xbar - mu0) / se
        if alternative == "two-sided":
            # 작은 쪽 꼬리를 골라 두 배 한다. z의 부호를 따지지 않아도 되고
            # 어느 쪽으로 치우쳐도 같은 식이 쓰인다.
            p = 2 * min(norm.cdf(z), norm.sf(z))
        elif alternative == "less":
            p = norm.cdf(z)
        else:
            # 오른쪽 꼬리는 sf로 계산한다. 1 - cdf는 꼬리에서 정밀도를 잃는다.
            p = norm.sf(z)
        return z, p, (p < alpha), "z-test"

    if sd is None:
        raise ValueError("Provide sd for t-test or known_sigma for z-test.")
    se = sd / math.sqrt(n)
    df = n - 1               # sd를 자료에서 추정했으므로 자유도 하나를 잃는다
    t = (xbar - mu0) / se
    if alternative == "two-sided":
        p = 2 * min(tdist.cdf(t, df), tdist.sf(t, df))
    elif alternative == "less":
        p = tdist.cdf(t, df)
    else:
        p = tdist.sf(t, df)
    return t, p, (p < alpha), f"t-test (df={df})"

import numpy as np
from scipy import stats

rng = np.random.default_rng(9)
x = rng.normal(3.4, 1.2, 25)
n, xbar, sd = len(x), x.mean(), x.std(ddof=1)
mu0 = 3.0

print(f"n = {n},  xbar = {xbar:.4f},  s = {sd:.4f}")
for alt in ("two-sided", "greater", "less"):
    t_my, p_my, _, _ = test_mean_one_sample(xbar, n, mu0, sd=sd, alternative=alt)
    r = stats.ttest_1samp(x, mu0, alternative=alt)
    print(f"  {alt:>9s}:  요약통계량 p = {p_my:.8f}   scipy p = {r.pvalue:.8f}")

# 대칭 항등식 2 min(F(t), 1-F(t)) = P(|T| >= |t|) 를 여러 t 에서 확인한다.
grid = np.array([-3.1, -0.4, 0.0, 0.9, 2.7])
tail = 2 * stats.t.sf(np.abs(grid), n - 1)                                 # P(|T| >= |t|)
mins = 2 * np.minimum(stats.t.cdf(grid, n - 1), stats.t.sf(grid, n - 1))   # 함수가 쓰는 식
print(f"두 식의 최대 차이 = {np.max(np.abs(tail - mins)):.2e}")

n = 25,  xbar = 3.5836,  s = 1.3921
  two-sided:  요약통계량 p = 0.04680592   scipy p = 0.04680592
    greater:  요약통계량 p = 0.02340296   scipy p = 0.02340296
       less:  요약통계량 p = 0.97659704   scipy p = 0.97659704
두 식의 최대 차이 = 0.00e+00

stat, p, reject, label = test_mean_one_sample(
    xbar=3.2, n=25, mu0=3.0, sd=1.1, alternative="greater"
)
print(label, "stat:", stat, "p:", p, "reject:", reject)

# 같은 자료를 sigma=1.1을 안다고 가정하고 z-검정으로도 해 본다.
stat_z, p_z, reject_z, label_z = test_mean_one_sample(
    xbar=3.2, n=25, mu0=3.0, known_sigma=1.1, alternative="greater"
)
print(label_z, "stat:", stat_z, "p:", p_z, "reject:", reject_z)

t-test (df=24) stat: 0.9090909090909097 p: 0.18617076763866552 reject: False
z-test stat: 0.9090909090909097 p: 0.1816510704434488 reject: False

import math

import numpy as np
from scipy import stats
from scipy.special import gammaln

nu = 24
t_obs = 0.2 / (1.1 / math.sqrt(25))

# c_nu = E[sqrt(W/nu)],  W ~ chi2_nu.  1 보다 작다는 것이 요점이다.
c_nu = math.exp(0.5 * math.log(2 / nu) + gammaln((nu + 1) / 2) - gammaln(nu / 2))

print(f"t_obs = {t_obs:.6f},  c_24 = {c_nu:.6f}")
print(f"정규 꼬리   Phibar(t)      = {stats.norm.sf(t_obs):.6f}")
print(f"옌센 하한   Phibar(c t)    = {stats.norm.sf(c_nu * t_obs):.6f}")
print(f"t 꼬리      P(T_24 >= t)   = {stats.t.sf(t_obs, nu):.6f}")

# 부등식이 t 전체에서 성립하는지 격자로 훑는다.
tg = np.linspace(0.01, 8.0, 200_001)
gap = stats.t.sf(tg, nu) - stats.norm.sf(tg)
print(f"(t 꼬리 - 정규 꼬리) 의 최솟값 = {gap.min():.3e}")

# sigma 를 안다고 잘못 가정하면 실제 유의수준이 얼마가 되는가.
z975 = stats.norm.ppf(0.975)
print("\n    n      nu   실제 크기   올바른 임계값")
for n in (2, 3, 5, 11, 25, 61, 101, 1001):
    print(f"{n:5d}  {n - 1:6d}      {2 * stats.t.sf(z975, n - 1):.4f}"
          f"          {stats.t.ppf(0.975, n - 1):.3f}")

t_obs = 0.909091,  c_24 = 0.989640
정규 꼬리   Phibar(t)      = 0.181651
옌센 하한   Phibar(c t)    = 0.184147
t 꼬리      P(T_24 >= t)   = 0.186171
(t 꼬리 - 정규 꼬리) 의 최솟값 = 1.578e-08

    n      nu   실제 크기   올바른 임계값
    2       1      0.3003          12.706
    3       2      0.1891          4.303
    5       4      0.1216          2.776
   11      10      0.0784          2.228
   25      24      0.0617          2.064
   61      60      0.0546          2.000
  101     100      0.0528          1.984
 1001    1000      0.0503          1.962

import numpy as np
from scipy import stats

rng = np.random.default_rng(19)
x = np.array([12.4, 9.8, 14.1, 11.2, 13.7, 10.5, 15.3,
              12.9, 11.8, 13.1, 10.2, 14.6])
mu0 = 11.0
n, B = len(x), 99_999

# ① t 검정
t_obs = (x.mean() - mu0) / (x.std(ddof=1) / np.sqrt(n))
print(f"t({n-1}) = {t_obs:.4f},  p = {2 * stats.t.sf(abs(t_obs), n-1):.4f}")

# ② 부트스트랩-t (자료를 mu0 로 중심 이동한 뒤 재표본)
xc = x - x.mean() + mu0
idx = rng.integers(0, n, (B, n))
xs = xc[idx]
tb = (xs.mean(1) - mu0) / (xs.std(1, ddof=1) / np.sqrt(n))
print(f"부트스트랩 p = {(np.sum(np.abs(tb) >= abs(t_obs)) + 1) / (B + 1):.4f}")

# ③ 부호 뒤집기 순열
d = x - mu0
s = rng.choice([-1.0, 1.0], (B, n))
ts = (d * s).mean(1) / ((d * s).std(1, ddof=1) / np.sqrt(n))
print(f"부호 순열 p  = {(np.sum(np.abs(ts) >= abs(t_obs)) + 1) / (B + 1):.4f}")

# ④ 윌콕슨
print(f"윌콕슨 p     = {stats.wilcoxon(d).pvalue:.4f}")

t(11) = 2.8271,  p = 0.0165
부트스트랩 p = 0.0202
부호 순열 p  = 0.0192
윌콕슨 p     = 0.0269

rng = np.random.default_rng(21)
M, n = 4_000, 20
cnt_t = cnt_perm = 0
for _ in range(M):
    y = rng.exponential(1.0, n)          # 평균 1, 오른쪽으로 치우침
    t = (y.mean() - 1.0) / (y.std(ddof=1) / np.sqrt(n))
    cnt_t += abs(t) > stats.t.ppf(0.975, n - 1)
    d = y - 1.0
    s = rng.choice([-1.0, 1.0], (999, n))
    ts = (d * s).mean(1) / ((d * s).std(1, ddof=1) / np.sqrt(n))
    cnt_perm += (np.sum(np.abs(ts) >= abs(t)) + 1) / 1000 < 0.05
print(f"지수분포(n=20): t 검정 {cnt_t / M:.4f},  "
      f"부호 순열 {cnt_perm / M:.4f}")

지수분포(n=20): t 검정 0.0862,  부호 순열 0.0860

import numpy as np
from scipy import stats

x = np.array([52.1, 48.7, 51.3, 49.8, 53.2, 50.4, 52.8,
              49.1, 51.7, 50.9, 48.3, 52.5])
mu0 = 50.0

def report(v, tag):
    n = len(v)
    t = stats.ttest_1samp(v, mu0)
    w = stats.wilcoxon(v - mu0)
    k = (v > mu0).sum()
    sg = stats.binomtest(k, n, 0.5).pvalue
    tm = stats.trim_mean(v, 0.1)
    print(f"{tag:14s} 평균 {v.mean():7.3f}  t p={t.pvalue:.4f}  "
          f"윌콕슨 p={w.pvalue:.4f}  부호 p={sg:.4f}  "
          f"10% 절사평균 {tm:.3f}")

report(x, "원자료")
y = x.copy()
y[0] = 20.0                       # 기록 오류 하나
report(y, "이상치 1개")

원자료            평균  50.900  t p=0.0857  윌콕슨 p=0.1099  부호 p=0.3877  10% 절사평균 50.930
이상치 1개         평균  48.225  t p=0.5101  윌콕슨 p=0.5186  부호 p=0.7744  10% 절사평균 50.550

import numpy as np
from scipy import stats

rng = np.random.default_rng(77)
M, n = 20_000, 30

cnt_ok = 0       # μ0 를 사전에 고정한 경우
cnt_bad = 0      # μ0 를 같은 자료에서 고른 경우
for _ in range(M):
    x = rng.normal(0, 1, n)
    cnt_ok += stats.ttest_1samp(x, 0.0).pvalue < 0.05
    mu0_bad = np.percentile(x, 25)          # 자료에서 고른 μ0
    cnt_bad += stats.ttest_1samp(x, mu0_bad).pvalue < 0.05
print(f"사전 고정 μ0: 제1종 오류율 {cnt_ok / M:.4f}")
print(f"자료에서 고른 μ0: 기각률 {cnt_bad / M:.4f}")

사전 고정 μ0: 제1종 오류율 0.0511
자료에서 고른 μ0: 기각률 0.9769

import numpy as np
from scipy import stats

rng = np.random.default_rng(31)
M = 40_000
print(f"{'분포':>8s} " + " ".join(f"{'n='+str(n):>8s}"
                                  for n in [10, 100, 1000]))
for name, gen, ok in [("정규", lambda s: rng.normal(0, 1, s), True),
                      ("t(2.5)", lambda s: rng.standard_t(2.5, s), True),
                      ("코시", lambda s: rng.standard_cauchy(s), False)]:
    row = []
    for n in [10, 100, 1000]:
        x = gen((M, n))
        t = x.mean(1) / (x.std(1, ddof=1) / np.sqrt(n))
        row.append(np.mean(np.abs(t) > 1.96))
    print(f"{name:>8s} " + " ".join(f"{v:8.4f}" for v in row))

      분포     n=10    n=100   n=1000
      정규   0.0838   0.0512   0.0512
  t(2.5)   0.0695   0.0477   0.0490
      코시   0.0407   0.0217   0.0209