일표본 비율 검정¶

from scipy.stats import norm, binomtest
import math

def test_prop_one_sample(k, n, p0=0.5, method="score",
                         alt="two-sided", alpha=0.05):
    """method='score'(정규근사) 또는 'exact'(이항).

    정확검정에는 검정통계량이 따로 없다. 이항분포에서 꼬리확률을
    바로 더하므로 표준화할 대상이 없기 때문이다. 그래서 None을 돌려준다.
    """
    phat = k / n
    if method == "exact":
        p = binomtest(k, n, p0, alternative=alt).pvalue
        return None, p, (p < alpha), "exact binomial"

    # 표준오차에 phat이 아니라 **p0**을 넣는다.
    # H0가 참이라는 가정 아래의 확률을 재는 것이 검정이기 때문이다.
    # 신뢰구간에서는 phat을 넣었다. 목적이 다르면 대입하는 값도 다르다.
    se0 = math.sqrt(p0 * (1 - p0) / n)
    z = (phat - p0) / se0
    if alt == "two-sided":
        p = 2 * min(norm.cdf(z), 1 - norm.cdf(z))
    elif alt == "less":
        p = norm.cdf(z)
    else:
        p = 1 - norm.cdf(z)
    return z, p, (p < alpha), "score z-test"

from fractions import Fraction

z975 = norm.ppf(0.975)


def accept_interval(n, p0, z=z975):
    """점수 검정이 채택하는 k 의 구간 [np0 - c, np0 + c]."""
    c = z * math.sqrt(n * p0 * (1 - p0))
    return n * p0 - c, n * p0 + c


def exact_size(reject, n, p0):
    """기각역의 확률을 유리수로 정확히 더한다. 부동소수점 오차가 없다."""
    q = Fraction(p0).limit_denominator()
    return sum(Fraction(math.comb(n, k)) * q**k * (1 - q)**(n - k)
               for k in reject)


n, p0 = 5, 0.5
lo, hi = accept_interval(n, p0)
rej_score = [k for k in range(n + 1) if not (lo <= k <= hi)]
rej_exact = [k for k in range(n + 1) if binomtest(k, n, p0).pvalue <= 0.05]

print(f"n = {n}, p0 = {p0}:  채택 구간 = [{lo:.6f}, {hi:.6f}]")
print(f"  점수 검정 기각역 = {rej_score},  정확 크기 = "
      f"{exact_size(rej_score, n, p0)} = {float(exact_size(rej_score, n, p0)):.4f}")
print(f"  정확검정 기각역 = {rej_exact},  정확 크기 = "
      f"{exact_size(rej_exact, n, p0)}")
print("  k 마다의 정확 양측 p-값:",
      [float(binomtest(k, n, p0).pvalue) for k in range(n + 1)])

print("\np0 = 0.5, 명목 0.05 에서 두 검정의 정확한 크기")
print("   n  점수 기각역                   크기     정확검정 기각역        크기")
for n in range(5, 15):
    lo, hi = accept_interval(n, p0)
    rs = [k for k in range(n + 1) if not (lo <= k <= hi)]
    re_ = [k for k in range(n + 1) if binomtest(k, n, p0).pvalue <= 0.05]
    print(f"{n:4d}  {str(rs):28s} {float(exact_size(rs, n, p0)):.4f}   "
          f"{str(re_):20s} {float(exact_size(re_, n, p0)):.4f}")

n = 5, p0 = 0.5:  채택 구간 = [0.308694, 4.691306]
  점수 검정 기각역 = [0, 5],  정확 크기 = 1/16 = 0.0625
  정확검정 기각역 = [],  정확 크기 = 0
  k 마다의 정확 양측 p-값: [0.0625, 0.375, 1.0, 1.0, 0.375, 0.0625]

p0 = 0.5, 명목 0.05 에서 두 검정의 정확한 크기
   n  점수 기각역                   크기     정확검정 기각역        크기
   5  [0, 5]                       0.0625   []                   0.0000
   6  [0, 6]                       0.0312   [0, 6]               0.0312
   7  [0, 7]                       0.0156   [0, 7]               0.0156
   8  [0, 1, 7, 8]                 0.0703   [0, 8]               0.0078
   9  [0, 1, 8, 9]                 0.0391   [0, 1, 8, 9]         0.0391
  10  [0, 1, 9, 10]                0.0215   [0, 1, 9, 10]        0.0215
  11  [0, 1, 2, 9, 10, 11]         0.0654   [0, 1, 10, 11]       0.0117
  12  [0, 1, 2, 10, 11, 12]        0.0386   [0, 1, 2, 10, 11, 12] 0.0386
  13  [0, 1, 2, 11, 12, 13]        0.0225   [0, 1, 2, 11, 12, 13] 0.0225
  14  [0, 1, 2, 3, 11, 12, 13, 14] 0.0574   [0, 1, 2, 12, 13, 14] 0.0129

stat, p, reject, label = test_prop_one_sample(
    k=12, n=50, p0=0.2, method="score", alt="two-sided"
)
print(label, "stat:", stat, "p:", p, "reject:", reject)

# 같은 자료의 정확한 이항검정
stat_e, p_e, reject_e, label_e = test_prop_one_sample(
    k=12, n=50, p0=0.2, method="exact", alt="two-sided"
)
print(label_e, "stat:", stat_e, "p:", p_e, "reject:", reject_e)

score z-test stat: 0.7071067811865471 p: 0.4795001221869537 reject: False
exact binomial stat: None p: 0.47974220659401984 reject: False

import math
from fractions import Fraction

from scipy.stats import norm, binomtest


def four_pvalues(k, n, p0):
    """같은 자료에 대한 네 가지 양측 p-값."""
    q = Fraction(p0).limit_denominator()
    w = [Fraction(math.comb(n, j)) * q**j * (1 - q)**(n - j) for j in range(n + 1)]

    sd = math.sqrt(n * p0 * (1 - p0))
    z = (k - n * p0) / sd                      # 점수 z (k 단위로 적은 것)
    p_score = 2 * norm.sf(abs(z))
    z_cc = (abs(k - n * p0) - 0.5) / sd        # 경계를 반 칸 밖으로 민다
    p_cc = 2 * norm.sf(z_cc)

    # 정확 ①: 관측된 칸보다 확률이 작거나 같은 칸을 모두 더한다(scipy 의 기본).
    p_small = sum(w[j] for j in range(n + 1) if w[j] <= w[k])
    # 정확 ②: 작은 쪽 꼬리를 두 배 한다.
    p_twice = 2 * min(sum(w[: k + 1]), sum(w[k:]))
    return z, p_score, z_cc, p_cc, p_small, p_twice


for k, n, p0 in [(2, 10, 0.5), (12, 50, 0.2)]:
    z, p_score, z_cc, p_cc, p_small, p_twice = four_pvalues(k, n, p0)
    print(f"k = {k}, n = {n}, p0 = {p0}")
    print(f"  점수            z = {z:+.6f}  p = {p_score:.6f}")
    print(f"  연속성 보정    z = {z_cc:+.6f}  p = {p_cc:.6f}")
    frac = f"  = {p_small}" if p_small.denominator < 10**6 else ""
    print(f"  정확 (작은 확률 합)       p = {float(p_small):.6f}{frac}")
    print(f"  정확 (작은 꼬리 두 배)    p = {float(p_twice):.6f}")
    print(f"  scipy binomtest           p = {binomtest(k, n, p0).pvalue:.6f}")
    print(f"  두 정확 관례의 차이 = {float(p_twice - p_small):.6f}")

k = 2, n = 10, p0 = 0.5
  점수            z = -1.897367  p = 0.057780
  연속성 보정    z = +1.581139  p = 0.113846
  정확 (작은 확률 합)       p = 0.109375  = 7/64
  정확 (작은 꼬리 두 배)    p = 0.109375
  scipy binomtest           p = 0.109375
  두 정확 관례의 차이 = 0.000000
k = 12, n = 50, p0 = 0.2
  점수            z = +0.707107  p = 0.479500
  연속성 보정    z = +0.530330  p = 0.595883
  정확 (작은 확률 합)       p = 0.479742
  정확 (작은 꼬리 두 배)    p = 0.578665
  scipy binomtest           p = 0.479742
  두 정확 관례의 차이 = 0.098923

import numpy as np
from scipy import stats

def levels(n, p0, alpha=0.05):
    k = np.arange(n + 1)
    ph = k / n
    z = stats.norm.ppf(1 - alpha / 2)
    se_w = np.sqrt(np.where((ph > 0) & (ph < 1), ph * (1 - ph) / n, np.inf))
    wald = np.abs(ph - p0) / se_w > z
    score = np.abs(ph - p0) / np.sqrt(p0 * (1 - p0) / n) > z
    a = np.where(k > 0, k * np.log(np.where(k > 0, ph / p0, 1)), 0)
    b = np.where(k < n, (n - k) * np.log(np.where(k < n,
                                                  (1 - ph) / (1 - p0), 1)), 0)
    lr = 2 * (a + b) > stats.chi2.ppf(1 - alpha, 1)
    ex = np.array([stats.binomtest(int(j), n, p0).pvalue
                   for j in k]) <= alpha
    pmf = stats.binom.pmf(k, n, p0)
    return [pmf[m].sum() for m in (wald, score, lr, ex)]

print(f"{'n':>5s} {'p0':>6s} {'왈드':>8s} {'점수':>8s} "
      f"{'우도비':>8s} {'정확':>8s}")
for n, p0 in [(20, 0.5), (20, 0.1), (50, 0.5), (50, 0.1),
              (100, 0.05), (200, 0.02)]:
    print(f"{n:5d} {p0:6.2f} " + " ".join(f"{v:8.4f}"
                                          for v in levels(n, p0)))

    n     p0       왈드       점수      우도비       정확
   20   0.50   0.0414   0.0414   0.0414   0.0414
   20   0.10   0.0024   0.0432   0.1328   0.0432
   50   0.50   0.0649   0.0649   0.0649   0.0328
   50   0.10   0.1159   0.0297   0.0583   0.0297
  100   0.05   0.1166   0.0341   0.0653   0.0341
  200   0.02   0.0743   0.0669   0.0378   0.0378

import numpy as np
from scipy import stats

def exact_power(n, p0, p1, alpha=0.05):
    k = np.arange(n + 1)
    pv = np.array([stats.binomtest(int(j), n, p0).pvalue for j in k])
    rej = pv <= alpha
    return (stats.binom.pmf(k, n, p1)[rej].sum(),
            stats.binom.pmf(k, n, p0)[rej].sum())

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"{'p0':>6s} {'p1':>6s} {'정확 n':>8s} {'근사 n':>8s} "
      f"{'근사 n 의 실제 검정력':>22s}")
for p0, p1 in [(0.5, 0.7), (0.5, 0.6), (0.1, 0.25), (0.05, 0.15),
               (0.02, 0.05)]:
    n_app = int(np.ceil((za * np.sqrt(p0 * (1 - p0))
                         + zb * np.sqrt(p1 * (1 - p1)))**2
                        / (p1 - p0)**2))
    n_ex = next(m for m in range(5, 3000)
                if exact_power(m, p0, p1)[0] >= 0.80)
    print(f"{p0:6.2f} {p1:6.2f} {n_ex:8d} {n_app:8d} "
          f"{exact_power(n_app, p0, p1)[0]:22.4f}")

    p0     p1     정확 n     근사 n          근사 n 의 실제 검정력
  0.50   0.70       49       47                 0.7801
  0.50   0.60      199      194                 0.7643
  0.10   0.25       44       41                 0.7296
  0.05   0.15       52       53                 0.8262
  0.02   0.05      226      233                 0.7322

for n in range(45, 56):
    print(f"n={n}: 수준 {exact_power(n, 0.5, 0.7)[1]:.4f}  "
          f"검정력 {exact_power(n, 0.5, 0.7)[0]:.4f}")

n=45: 수준 0.0357  검정력 0.7462
n=46: 수준 0.0259  검정력 0.7128
n=47: 수준 0.0400  검정력 0.7801
n=48: 수준 0.0293  검정력 0.7495
n=49: 수준 0.0444  검정력 0.8100
n=50: 수준 0.0328  검정력 0.7822
n=51: 수준 0.0489  검정력 0.8363
n=52: 수준 0.0365  검정력 0.8112
n=53: 수준 0.0270  검정력 0.7844
n=54: 수준 0.0402  검정력 0.8368
n=55: 수준 0.0300  검정력 0.8125

import numpy as np
from scipy import stats

n, m, rho = 1_000, 20, 0.05          # 군집 50개 × 20명
k, p0 = 560, 0.5
ph = k / n
deff = 1 + (m - 1) * rho

for name, d in [("독립 가정", 1.0), ("DEFF 반영", deff)]:
    se = np.sqrt(p0 * (1 - p0) * d / n)
    z = (ph - p0) / se
    print(f"{name:10s} SE {se:.5f}  z {z:6.3f}  "
          f"p {2 * stats.norm.sf(abs(z)):.5f}")
print(f"\nDEFF = {deff:.2f},  유효표본 {n / deff:.0f}명")

독립 가정      SE 0.01581  z  3.795  p 0.00015
DEFF 반영    SE 0.02208  z  2.717  p 0.00658

DEFF = 1.95,  유효표본 513명

rng = np.random.default_rng(7)
for cv in [0.0, 0.3, 0.6, 1.0]:
    w = 1 + cv * rng.standard_normal(1_000)
    w = np.clip(w, 0.1, None)
    neff = w.sum()**2 / (w**2).sum()
    print(f"가중치 변동계수 {cv:.1f}: 유효표본 {neff:.0f}명 "
          f"(설계효과 {1000 / neff:.2f})")

가중치 변동계수 0.0: 유효표본 1000명 (설계효과 1.00)
가중치 변동계수 0.3: 유효표본 913명 (설계효과 1.10)
가중치 변동계수 0.6: 유효표본 757명 (설계효과 1.32)
가중치 변동계수 1.0: 유효표본 636명 (설계효과 1.57)

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

n, k, p0 = 100, 61, 0.5

# ① 빈도주의
print(f"정확 이항 p = {stats.binomtest(k, n, p0).pvalue:.4f}")
ph = k / n
z = (ph - p0) / np.sqrt(p0 * (1 - p0) / n)
print(f"점수 z = {z:.4f},  p = {2 * stats.norm.sf(abs(z)):.4f}")

# ② 사후확률 (제프리스 사전분포)
a, b = k + 0.5, n - k + 0.5
print(f"\n사후분포 Beta({a}, {b})")
print(f"P(p > 0.5 | 자료) = {stats.beta.sf(0.5, a, b):.4f}")
lo, hi = stats.beta.ppf([0.025, 0.975], a, b)
print(f"95% 신용구간 ({lo:.4f}, {hi:.4f})")

# ③ 베이즈 인자 (H0: p=0.5 대 H1: p ~ Beta(1,1))
bf = np.exp(betaln(k + 1, n - k + 1) - betaln(1, 1) - n * np.log(0.5))
print(f"\nBF10 = {bf:.4f}  →  사전 1:1 이면 H1 의 사후확률 "
      f"{bf / (1 + bf):.4f}")

# ④ ROPE 접근: |p - 0.5| < 0.05 를 '실질적으로 같음'으로
rope = stats.beta.cdf(0.55, a, b) - stats.beta.cdf(0.45, a, b)
print(f"P(|p - 0.5| < 0.05 | 자료) = {rope:.4f}")

정확 이항 p = 0.0352
점수 z = 2.2000,  p = 0.0278

사후분포 Beta(61.5, 39.5)
P(p > 0.5 | 자료) = 0.9863
95% 신용구간 (0.5124, 0.7014)

BF10 = 1.3924  →  사전 1:1 이면 H1 의 사후확률 0.5820

P(|p - 0.5| < 0.05 | 자료) = 0.1128

import numpy as np
from scipy import stats

rng = np.random.default_rng(88)
M, p, days, per_day = 3_000, 0.10, 30, 200

cnt = 0
for _ in range(M):
    ka = kb = na = nb = 0
    for d in range(days):
        ka += rng.binomial(per_day, p)
        kb += rng.binomial(per_day, p)     # 참 차이 없음
        na += per_day
        nb += per_day
        pa, pb = ka / na, kb / nb
        pbar = (ka + kb) / (na + nb)
        se = np.sqrt(pbar * (1 - pbar) * (1 / na + 1 / nb))
        if se > 0 and abs(pa - pb) / se > 1.96:
            cnt += 1
            break
print(f"매일 확인하며 유의하면 중단: 거짓양성률 {cnt / M:.4f}")

매일 확인하며 유의하면 중단: 거짓양성률 0.2730

alpha1 = 2 * stats.norm.sf(3)
print(f"3σ 한계의 개별 오류율 {alpha1:.6f}")
for n_pts in [20, 50, 100, 500]:
    print(f"  {n_pts:4d}점 관측 시 한 번이라도 경보: "
          f"{1 - (1 - alpha1)**n_pts:.4f}")

3σ 한계의 개별 오류율 0.002700
    20점 관측 시 한 번이라도 경보: 0.0526
    50점 관측 시 한 번이라도 경보: 0.1264
   100점 관측 시 한 번이라도 경보: 0.2369
   500점 관측 시 한 번이라도 경보: 0.7412