기각역 시연¶

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

np.random.seed(42)
data = stats.norm.rvs(loc=170, scale=8, size=250)

mu0 = 172
n = len(data)
df = n - 1
xbar = data.mean()
s = data.std(ddof=1)
se = s / np.sqrt(n)

alpha = 0.05
# 양측이므로 alpha를 두 꼬리에 반씩 나눈다. ppf에 1 - alpha/2를 넣는 이유다.
t_crit = stats.t.ppf(1 - alpha / 2, df)
t_stat = (xbar - mu0) / se

print(f"x-bar = {xbar:.2f}, SE = {se:.2f}")
print(f"t-stat = {t_stat:.4f}, t-crit = +/-{t_crit:.4f}")
# 기각역을 t 척도가 아니라 cm 척도로도 적어 둔다.
# 실무자에게는 "표본평균이 171.04 아래면 기각"이 t보다 읽기 쉽다.
print(f"Rejection boundaries: {mu0 - t_crit*se:.2f} and {mu0 + t_crit*se:.2f}")
print(f"더 정밀하게: x-bar = {xbar:.6f}, s = {s:.6f}, SE = {se:.6f}")
print(f"             t = {t_stat:.6f},  경계 "
      f"{mu0 - t_crit * se:.6f} / {mu0 + t_crit * se:.6f}")
print(f"             p = {2 * stats.t.sf(abs(t_stat), df):.3e}")

# (2) 반너비가 2 cm 아래로 내려가는 첫 n 을 찾는다.
# **격자에서 argmin/argmax 로 찾지 말 것.** 조건을 만족하는 칸이 하나도
# 없어도 0 을 돌려주므로, 조건을 만족하는 n 이 실제로 있는지부터 확인한다.
ns = np.arange(10, 400)
half = stats.t.ppf(0.975, ns - 1) * s / np.sqrt(ns)
ok = np.flatnonzero(half <= 2.0)
assert ok.size > 0, "반너비가 2 cm 이하가 되는 n 이 격자 안에 없다"
n2 = ns[ok[0]]
print(f"\n반너비 <= 2 cm 인 첫 n = {n2}  (반너비 {half[ok[0]]:.4f} cm)"
      f"   바로 앞 n = {n2 - 1} 은 {half[ok[0] - 1]:.4f} cm")
print(f"z 어림: n = {(stats.norm.ppf(0.975) * s / 2) ** 2:.1f}")

x-bar = 169.98, SE = 0.49
t-stat = -4.1314, t-crit = +/-1.9695
Rejection boundaries: 171.04 and 172.96
더 정밀하게: x-bar = 169.980617, s = 7.728430, SE = 0.488789
             t = -4.131402,  경계 171.037312 / 172.962688
             p = 4.927e-05

반너비 <= 2 cm 인 첫 n = 60  (반너비 1.9965 cm)   바로 앞 n = 59 은 2.0140 cm
z 어림: n = 57.4

# 기각역은 자료를 보기 전에 정해지는 영역이다. 관측된 통계량이 칠해진
# 구역 안에 떨어지면 기각한다. p-값 방식과 결론은 언제나 같다.
x_t = np.linspace(-5, 5, 300)
y_t = stats.t.pdf(x_t, df)

# 그림에서 읽을 수치를 미리 찍어 둔다.
print(f"꼭대기 밀도 f(0)        = {stats.t.pdf(0, df):.6f}")
print(f"임계값에서의 밀도 f(tc) = {stats.t.pdf(t_crit, df):.6f}")
print(f"관측값에서의 밀도 f(t)  = {stats.t.pdf(t_stat, df):.8f}"
      f"   꼭대기의 1/{stats.t.pdf(0, df) / stats.t.pdf(t_stat, df):.0f}")
print(f"칠한 꼬리 넓이 (한쪽)   = {stats.t.sf(t_crit, df):.6f}  (두쪽 합 "
      f"{2 * stats.t.sf(t_crit, df):.4f})")
print(f"그림 밖 |t| > 5 의 질량 = {2 * stats.t.sf(5, df):.3e}")
print(f"P(|T| > |t|) = {2 * stats.t.sf(abs(t_stat), df):.3e}"
      f"   칠한 왼쪽 꼬리의 {stats.t.cdf(t_stat, df) / stats.t.cdf(-t_crit, df):.4%}")

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(x_t, y_t, "tomato", lw=2)
ax.fill_between(x_t[x_t <= -t_crit], stats.t.pdf(x_t[x_t <= -t_crit], df),
                color="tomato", alpha=0.5, label="Rejection region")
ax.fill_between(x_t[x_t >= t_crit], stats.t.pdf(x_t[x_t >= t_crit], df),
                color="tomato", alpha=0.5)
ax.axvline(t_stat, color="blue", linestyle="--", label=f"t = {t_stat:.2f}")
ax.set_xlabel("t")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()

꼭대기 밀도 f(0)        = 0.398542
임계값에서의 밀도 f(tc) = 0.057714
관측값에서의 밀도 f(t)  = 0.00010028   꼭대기의 1/3974
칠한 꼬리 넓이 (한쪽)   = 0.025000  (두쪽 합 0.0500)
그림 밖 |t| > 5 의 질량 = 1.082e-06
P(|T| > |t|) = 4.927e-05   칠한 왼쪽 꼬리의 0.0985%

df = 100
alpha = 0.05
x = np.linspace(-5, 5, 300)
y = stats.t.pdf(x, df)

# 단측이므로 alpha를 나누지 않는다. 꼬리 하나에 5%를 통째로 준다.
t_lo = stats.t.ppf(alpha, df)
print(f"Left-tailed critical value: {t_lo:.4f}")

t_hi = stats.t.ppf(1 - alpha, df)
print(f"Right-tailed critical value: {t_hi:.4f}")

t_two = stats.t.ppf(1 - alpha / 2, df)
print(f"Two-tailed critical value : +/-{t_two:.4f}")
print(f"대칭성 확인: t_lo + t_hi = {t_lo + t_hi:.1e}")

def power_one(lam):      # 우측 단측검정
    return stats.nct.sf(t_hi, df, lam)

def power_two(lam):      # 양측검정
    return stats.nct.sf(t_two, df, lam) + stats.nct.cdf(-t_two, df, lam)

print(f"\n{'lambda':>7} {'우측 단측':>10} {'양측':>8} {'차이':>9}")
for lam in (-3.0, -1.0, 0.0, 1.0, 2.0, 3.0):
    p1, p2 = power_one(lam), power_two(lam)
    print(f"{lam:>7.1f} {p1:>10.4f} {p2:>8.4f} {p1 - p2:>+9.4f}")

Left-tailed critical value: -1.6602
Right-tailed critical value: 1.6602
Two-tailed critical value : +/-1.9840
대칭성 확인: t_lo + t_hi = -6.7e-16

 lambda      우측 단측       양측        차이
   -3.0     0.0000   0.8440   -0.8440
   -1.0     0.0042   0.1677   -0.1636
    0.0     0.0500   0.0500   -0.0000
    1.0     0.2573   0.1677   +0.0896
    2.0     0.6337   0.5083   +0.1253
    3.0     0.9090   0.8440   +0.0650

import numpy as np
from scipy import stats

for nu in [5, 10, 20, 50]:
    lo = stats.chi2.ppf(0.025, nu)
    hi = stats.chi2.ppf(0.975, nu)
    print(f"ν = {nu:3d}:  ({lo:7.3f}, {hi:7.3f})   "
          f"중심 {nu}에서의 거리 {nu - lo:6.3f} / {hi - nu:6.3f}   "
          f"밀도 {stats.chi2.pdf(lo, nu):.5f} / {stats.chi2.pdf(hi, nu):.5f}")

ν =   5:  (  0.831,  12.833)   거리  4.169 /  7.833   밀도 0.06651 / 0.00999
ν =  10:  (  3.247,  20.483)   거리  6.753 / 10.483   밀도 0.02854 / 0.00817
ν =  20:  (  9.591,  34.170)   거리 10.409 / 14.170   밀도 0.01528 / 0.00650
ν =  50:  ( 32.357,  71.420)   거리 17.643 / 21.420   밀도 0.00784 / 0.00462

import numpy as np
from scipy import stats

n, p0, alpha = 20, 0.5, 0.05
a_hi = stats.binom.sf(14, n, p0)          # P(X ≥ 15)
pmf14 = stats.binom.pmf(14, n, p0)
gamma = (alpha - a_hi) / pmf14
print(f"P(X ≥ 15) = {a_hi:.5f}")
print(f"P(X = 14)  = {pmf14:.5f}")
print(f"γ = {gamma:.4f}   →  수준 {a_hi + gamma * pmf14:.5f}")

# 검정력 비교
for p1 in [0.6, 0.7, 0.8]:
    pw0 = stats.binom.sf(14, n, p1)
    pw1 = pw0 + gamma * stats.binom.pmf(14, n, p1)
    print(f"p1 = {p1}: 비무작위화 {pw0:.4f}   무작위화 {pw1:.4f}")

P(X ≥ 15) = 0.02069
P(X = 14)  = 0.03696
γ = 0.7928   →  수준 0.05000
p1 = 0.6: 비무작위화 0.1256   무작위화 0.2242
p1 = 0.7: 비무작위화 0.4164   무작위화 0.5683
p1 = 0.8: 비무작위화 0.8042   무작위화 0.8907

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
M, p = 40_000, 2
z = stats.norm.ppf(1 - 0.05 / (2 * p))
crit = stats.chi2.ppf(0.95, p)

print(f"{'ρ':>5s} {'타원 수준':>10s} {'본페 수준':>10s} "
      f"{'타원 검정력':>12s} {'본페 검정력':>12s}")
for rho in [0.0, 0.5, 0.9]:
    S = np.array([[1, rho], [rho, 1]])
    L = np.linalg.cholesky(S)
    Si = np.linalg.inv(S)
    out = {}
    for d in [np.zeros(2), np.array([2.0, 0.0])]:
        x = rng.standard_normal((M, p)) @ L.T + d
        ell = np.einsum("ij,jk,ik->i", x, Si, x) > crit
        box = (np.abs(x) > z).any(1)
        out[d[0]] = (ell.mean(), box.mean())
    print(f"{rho:5.1f} {out[0.0][0]:10.4f} {out[0.0][1]:10.4f} "
          f"{out[2.0][0]:12.4f} {out[2.0][1]:12.4f}")

    ρ      타원 수준      본페 수준       타원 검정력       본페 검정력
  0.0     0.0497     0.0485       0.4197       0.4253
  0.5     0.0499     0.0464       0.5329       0.4175
  0.9     0.0511     0.0366       0.9895       0.4144

rng = np.random.default_rng(5)
print(f"{'p':>4s} {'타원 검정력':>12s} {'본페 검정력':>12s}")
for p in [2, 5, 10, 20]:
    z = stats.norm.ppf(1 - 0.05 / (2 * p))
    crit = stats.chi2.ppf(0.95, p)
    d = np.zeros(p)
    d[0] = 3.0                                # 한 좌표만 이동, 상관 없음
    x = rng.standard_normal((M, p)) + d
    print(f"{p:4d} {np.mean((x**2).sum(1) > crit):12.4f} "
          f"{np.mean((np.abs(x) > z).any(1)):12.4f}")

   p    타원 검정력    본페 검정력
   2       0.7703       0.7808
   5       0.6281       0.6799
  10       0.4931       0.5925
  20       0.3560       0.5144

import numpy as np
from scipy import stats

n = 10
nu = n - 1
lo = stats.chi2.ppf(0.025, nu)
hi = stats.chi2.ppf(0.975, nu)
print(f"기각역: W < {lo:.3f} 또는 W > {hi:.3f}")
print(f"{'σ²/σ0²':>8s} {'검정력':>10s}")
for r in [0.5, 0.7, 0.9, 0.95, 1.0, 1.05, 1.1, 1.3, 2.0]:
    pw = stats.chi2.cdf(lo / r, nu) + stats.chi2.sf(hi / r, nu)
    mark = "  ← α 미만" if pw < 0.05 else ""
    print(f"{r:8.2f} {pw:10.4f}{mark}")

기각역: W < 2.700 또는 W > 19.023
  σ²/σ0²        검정력
    0.50     0.2020
    0.70     0.0808
    0.90     0.0478  ← α 미만
    0.95     0.0476  ← α 미만
    1.00     0.0500
    1.05     0.0550
    1.10     0.0622
    1.30     0.1113
    2.00     0.3934

import numpy as np
from scipy import stats

n, s, alpha = 25, 15.0, 0.05
t = stats.t.ppf(1 - alpha / 2, n - 1)
se = s / np.sqrt(n)

print("가로 단면 — μ0 를 고정했을 때의 수용역 (x̄ 의 범위)")
for mu0 in [95, 100, 105]:
    print(f"  μ0 = {mu0:3d}:  ({mu0 - t * se:.3f}, {mu0 + t * se:.3f})")

print("\n세로 단면 — x̄ 를 고정했을 때의 신뢰구간 (μ 의 범위)")
for xbar in [98, 102, 107]:
    print(f"  x̄ = {xbar:3d}:  ({xbar - t * se:.3f}, {xbar + t * se:.3f})")

print("\n쌍대성 확인")
for mu0, xbar in [(100, 105), (100, 107), (105, 107)]:
    in_ci = abs(xbar - mu0) <= t * se
    print(f"  μ0={mu0}, x̄={xbar}:  구간에 {'포함' if in_ci else '미포함'}, "
          f"검정은 {'비기각' if in_ci else '기각'}")

가로 단면 — μ0 를 고정했을 때의 수용역 (x̄ 의 범위)
  μ0 =  95:  (88.808, 101.192)
  μ0 = 100:  (93.808, 106.192)
  μ0 = 105:  (98.808, 111.192)

세로 단면 — x̄ 를 고정했을 때의 신뢰구간 (μ 의 범위)
  x̄ =  98:  (91.808, 104.192)
  x̄ = 102:  (95.808, 108.192)
  x̄ = 107:  (100.808, 113.192)

쌍대성 확인
  μ0=100, x̄=105:  구간에 포함, 검정은 비기각
  μ0=100, x̄=107:  구간에 미포함, 검정은 기각
  μ0=105, x̄=107:  구간에 포함, 검정은 비기각