평균 차이에 대한 이표본 Z-검정¶

import numpy as np
from scipy import stats

def two_sample_z_test(x_bar, y_bar, sigma1, sigma2, n1, n2,
                      delta0=0, alternative="two-sided"):
    """두 평균의 차이에 대한 이표본 z 검정.

    Parameters
    ----------
    x_bar, y_bar : float
        Sample means.
    sigma1, sigma2 : float
        Known population standard deviations.
    n1, n2 : int
        Sample sizes.
    delta0 : float
        Hypothesized difference (default 0).
    alternative : str
        'two-sided', 'less', or 'greater'.

    Returns
    -------
    z_stat : float
        Test statistic.
    p_value : float
        P-value.
    """
    # sigma를 알고 있으므로 표본에서 추정하지 않는다.
    # 그래서 자유도라는 개념이 없고 표준정규를 그대로 쓴다.
    se = np.sqrt(sigma1**2 / n1 + sigma2**2 / n2)
    z_stat = ((x_bar - y_bar) - delta0) / se

    if alternative == "two-sided":
        p_value = 2 * (1 - stats.norm.cdf(abs(z_stat)))
    elif alternative == "less":
        p_value = stats.norm.cdf(z_stat)
    elif alternative == "greater":
        p_value = 1 - stats.norm.cdf(z_stat)
    else:
        raise ValueError("alternative must be 'two-sided', 'less', or 'greater'")

    return z_stat, p_value


# 보기 1 의 자료
z, p = two_sample_z_test(78, 74, 12, 15, 50, 60, alternative="two-sided")
print(f"z = {z:.3f}, p-value = {p:.3f}")

z = 1.553, p-value = 0.120

z95 = stats.norm.ppf(0.975)
sigma1, sigma2, n1, n2, diff = 12.0, 15.0, 50, 60, 4.0

# (1) 두 표본을 같은 배수 k 로 늘린다.
print("   k    n1    n2       SE        z         p")
for k in (1, 2, 3, 4):
    se = np.sqrt(sigma1**2 / (n1 * k) + sigma2**2 / (n2 * k))
    z = diff / se
    print(f"{k:4d} {n1 * k:5d} {n2 * k:5d} {se:8.4f} {z:8.4f}"
          f"  {2 * stats.norm.sf(z):.6f}")
z1 = diff / np.sqrt(sigma1**2 / n1 + sigma2**2 / n2)
print(f"\np < 0.05 가 되는 k 의 문턱 = (z_0.025/z_1)^2 = {(z95 / z1)**2:.6f}"
      f"  →  k = 2")

# (2) B 학군만 늘린다. n1 은 50 에 묶여 있다.
print("\n      n2       SE        z         p")
for m in (60, 176, 1000, 10**6):
    se = np.sqrt(sigma1**2 / n1 + sigma2**2 / m)
    z = diff / se
    print(f"{m:8d} {se:8.4f} {z:8.4f}  {2 * stats.norm.sf(z):.6f}")
se_inf = sigma1 / np.sqrt(n1)
print(f"{'무한':>8} {se_inf:8.4f} {diff / se_inf:8.4f}  "
      f"{2 * stats.norm.sf(diff / se_inf):.6f}   ← 하한")

need = sigma2**2 / ((diff / z95)**2 - sigma1**2 / n1)
print(f"\np < 0.05 에 필요한 n2 > {need:.4f}  →  n2 = {int(np.ceil(need))}")

   k    n1    n2       SE        z         p
   1    50    60   2.5749   1.5535  0.120311
   2   100   120   1.8207   2.1969  0.028025
   3   150   180   1.4866   2.6907  0.007130
   4   200   240   1.2874   3.1069  0.001890

p < 0.05 가 되는 k 의 문턱 = (z_0.025/z_1)^2 = 1.591804  →  k = 2

      n2       SE        z         p
      60   2.5749   1.5535  0.120311
     176   2.0392   1.9615  0.049816
    1000   1.7621   2.2700  0.023207
 1000000   1.6971   2.3569  0.018427
      무한   1.6971   2.3570  0.018422   ← 하한

p < 0.05 에 필요한 n2 > 175.0858  →  n2 = 176

import numpy as np
from scipy import stats

x1, s1, n1 = 82.0, 10.0, 50
x2, s2, n2 = 78.0, 12.0, 60
se = np.sqrt(s1**2 / n1 + s2**2 / n2)
z = (x1 - x2) / se
print(f"관측: SE={se:.4f}  z={z:.4f}  p={2 * stats.norm.sf(abs(z)):.4f}")

def power_z(n1, n2, delta, s1, s2, alpha=0.05):
    se = np.sqrt(s1**2 / n1 + s2**2 / n2)
    za = stats.norm.ppf(1 - alpha / 2)
    return (stats.norm.sf(za - delta / se)
            + stats.norm.cdf(-za - delta / se))

def n_z(delta, s1, s2, power=0.80, alpha=0.05, ratio=1.0):
    """n2 = ratio * n1 로 두었을 때 필요한 (n1, n2)."""
    za, zb = stats.norm.ppf(1 - alpha / 2), stats.norm.ppf(power)
    m = ((za + zb) / delta)**2 * (s1**2 + s2**2 / ratio)
    return int(np.ceil(m)), int(np.ceil(ratio * np.ceil(m)))

print(f"\n{'참 차이':>7s} {'현 설계 검정력':>14s} {'필요 n1':>8s} {'필요 n2':>8s}")
for d in [2, 3, 4, 5, 6]:
    a, b = n_z(d, s1, s2, ratio=1.2)
    print(f"{d:7.1f} {power_z(n1, n2, d, s1, s2):14.4f} {a:8d} {b:8d}")

print(f"\n총 {n1 + n2}명을 어떻게 나눌까 (참 차이 4)")
for a in [40, 45, 50, 55, 60]:
    b = n1 + n2 - a
    print(f"  n1={a:3d}, n2={b:3d}:  SE={np.sqrt(s1**2 / a + s2**2 / b):.4f}  "
          f"검정력 {power_z(a, b, 4.0, s1, s2):.4f}")
tot = n1 + n2
print(f"  네이만 최적 배분 σ1:σ2 = {s1:.0f}:{s2:.0f} → "
      f"n1={tot * s1 / (s1 + s2):.1f}, n2={tot * s2 / (s1 + s2):.1f}")

관측: SE=2.0976  z=1.9069  p=0.0565

   참 차이       현 설계 검정력    필요 n1    필요 n2
    2.0         0.1589      432      519
    3.0         0.2985      192      231
    4.0         0.4789      108      130
    5.0         0.6641       70       84
    6.0         0.8161       48       58

총 110명을 어떻게 나눌까 (참 차이 4)
  n1= 40, n2= 70:  SE=2.1347  검정력 0.4657
  n1= 45, n2= 65:  SE=2.1066  검정력 0.4757
  n1= 50, n2= 60:  SE=2.0976  검정력 0.4789
  n1= 55, n2= 55:  SE=2.1063  검정력 0.4758
  n1= 60, n2= 50:  SE=2.1323  검정력 0.4666
  네이만 최적 배분 σ1:σ2 = 10:12 → n1=50.0, n2=60.0

import numpy as np
from scipy import stats

x1, s1, n1 = 82.0, 10.0, 50
x2, s2, n2 = 78.0, 12.0, 60
se = np.sqrt(s1**2 / n1 + s2**2 / n2)
z = (x1 - x2) / se
print(f"z = {z:.4f}")
print(f"  양측        p = {2 * stats.norm.sf(abs(z)):.4f}  "
      f"→ {'기각' if 2 * stats.norm.sf(abs(z)) < 0.05 else '기각 못 함'}")
print(f"  단측 (μ1>μ2) p = {stats.norm.sf(z):.4f}  "
      f"→ {'기각' if stats.norm.sf(z) < 0.05 else '기각 못 함'}")
print(f"  단측 (μ1<μ2) p = {stats.norm.cdf(z):.4f}  "
      f"→ {'기각' if stats.norm.cdf(z) < 0.05 else '기각 못 함'}")

rng = np.random.default_rng(64)
M = 200_000
zz = rng.standard_normal(M)                 # H0 가 참일 때의 z
c = stats.norm.ppf(0.95)
print(f"\n사전에 방향을 지정한 단측검정: 제1종 오류율 {np.mean(zz > c):.4f}")
print(f"자료를 보고 방향을 고른 경우:   제1종 오류율 {np.mean(np.abs(zz) > c):.4f}")

z = 1.9069
  양측        p = 0.0565  → 기각 못 함
  단측 (μ1>μ2) p = 0.0283  → 기각
  단측 (μ1<μ2) p = 0.9717  → 기각 못 함

사전에 방향을 지정한 단측검정: 제1종 오류율 0.0506
자료를 보고 방향을 고른 경우:   제1종 오류율 0.1001

import numpy as np
from scipy import stats

x1, s1, n1 = 82.0, 10.0, 50
x2, s2, n2 = 78.0, 12.0, 60
diff = x1 - x2
se = np.sqrt(s1**2 / n1 + s2**2 / n2)
z95 = stats.norm.ppf(0.975)
lo, hi = diff - z95 * se, diff + z95 * se
print(f"차이 {diff:.1f},  SE {se:.4f},  95% CI ({lo:.4f}, {hi:.4f})\n")

for d0 in [0, 2, 4, 6, 8, 9]:
    zs = (diff - d0) / se
    p = 2 * stats.norm.sf(abs(zs))
    inside = lo <= d0 <= hi
    print(f"δ0={d0:2d}:  z={zs:7.4f}  p={p:.4f}  "
          f"{'기각 못 함' if p >= 0.05 else '기각':10s} "
          f"구간 안? {inside}")

print(f"\n비열등성 (Δ=2, '집단 1이 집단 2보다 2 이상 나쁘지는 않다')")
z_ni = (diff + 2) / se
print(f"  z={z_ni:.4f},  단측 p={stats.norm.sf(z_ni):.6f}  → 비열등 인정")

z90 = stats.norm.ppf(0.95)
lo90, hi90 = diff - z90 * se, diff + z90 * se
for D in [6.0, 8.0]:
    p_lo = stats.norm.sf((diff + D) / se)
    p_hi = stats.norm.cdf((diff - D) / se)
    tost = max(p_lo, p_hi)
    print(f"\n동등성 TOST (Δ={D:.0f}):  p = {tost:.4f}  "
          f"→ {'동등 인정' if tost < 0.05 else '동등 미인정'}")
    print(f"  90% CI ({lo90:.4f}, {hi90:.4f}) ⊂ (-{D:.0f}, {D:.0f})?  "
          f"{-D < lo90 and hi90 < D}")

차이 4.0,  SE 2.0976,  95% CI (-0.1113, 8.1113)

δ0= 0:  z= 1.9069  p=0.0565  기각 못 함     구간 안? True
δ0= 2:  z= 0.9535  p=0.3404  기각 못 함     구간 안? True
δ0= 4:  z= 0.0000  p=1.0000  기각 못 함     구간 안? True
δ0= 6:  z=-0.9535  p=0.3404  기각 못 함     구간 안? True
δ0= 8:  z=-1.9069  p=0.0565  기각 못 함     구간 안? True
δ0= 9:  z=-2.3837  p=0.0171  기각         구간 안? False

비열등성 (Δ=2, '집단 1이 집단 2보다 2 이상 나쁘지는 않다')
  z=2.8604,  단측 p=0.002116  → 비열등 인정

동등성 TOST (Δ=6):  p = 0.1702  → 동등 미인정
  90% CI (0.5497, 7.4503) ⊂ (-6, 6)?  False

동등성 TOST (Δ=8):  p = 0.0283  → 동등 인정
  90% CI (0.5497, 7.4503) ⊂ (-8, 8)?  True

import numpy as np
from scipy import stats

z95 = stats.norm.ppf(0.975)
print(f"{'가정 σ / 참 σ':>14s} {'구간 포함률':>12s} {'실제 α':>10s}")
for r in [0.6, 0.8, 0.9, 1.0, 1.1, 1.25, 1.5, 2.0]:
    cov = 2 * stats.norm.cdf(z95 * r) - 1
    print(f"{r:14.2f} {cov:12.4f} {1 - cov:10.4f}")

    가정 σ / 참 σ       구간 포함률       실제 α
          0.60       0.7604     0.2396
          0.80       0.8831     0.1169
          0.90       0.9223     0.0777
          1.00       0.9500     0.0500
          1.10       0.9689     0.0311
          1.25       0.9857     0.0143
          1.50       0.9967     0.0033
          2.00       0.9999     0.0001

from scipy import stats

z95 = stats.norm.ppf(0.975)
print(f"{'자유도':>7s} {'t 임계값':>10s} {'z 임계값':>10s} "
      f"{'z 를 쓸 때 실제 수준':>21s}")
for df in [5, 10, 20, 30, 50, 100, 200, 1000]:
    print(f"{df:7d} {stats.t.ppf(0.975, df):10.4f} {z95:10.4f} "
          f"{2 * stats.t.sf(z95, df):21.4f}")

    자유도      t 임계값      z 임계값         z 를 쓸 때 실제 수준
      5     2.5706     1.9600                0.1073
     10     2.2281     1.9600                0.0784
     20     2.0860     1.9600                0.0641
     30     2.0423     1.9600                0.0593
     50     2.0086     1.9600                0.0556
    100     1.9840     1.9600                0.0528
    200     1.9719     1.9600                0.0514
   1000     1.9623     1.9600                0.0503