베이즈 정리¶

import numpy as np

def bayes_theorem(prior, likelihood, evidence):
    """베이즈 정리:  사후 = (가능도 x 사전) / 증거

    P(H|E) = P(E|H) P(H) / P(E)
    """
    posterior = (likelihood * prior) / evidence
    return posterior

# 의학 진단 보기
prior_disease = 0.01         # 사전확률 P(질병) = 유병률 1%
sensitivity = 0.95           # 가능도 P(양성|질병) = 95%
specificity = 0.90           # 특이도 P(음성|건강) = 90%
false_positive_rate = 1 - specificity     # P(양성|건강) = 10%

# 분모(증거) P(양성)은 전확률의 법칙으로 구한다.
# 두 경로 — 질병이면서 양성, 건강하면서 양성 — 를 더한다.
p_positive = sensitivity * prior_disease + false_positive_rate * (1 - prior_disease)

# 이 분모가 베이즈 정리의 핵심이다. 건강한 사람이 99%나 되므로
# 위양성 10%가 만들어 내는 양성자 수가 진짜 환자보다 훨씬 많아진다.
posterior = bayes_theorem(prior_disease, sensitivity, p_positive)

print(f"P(disease | positive) = {posterior:.4f}")
print(f"Despite a 95% sensitive test, only {posterior*100:.1f}% of positives truly have the disease.")

# (1) 오즈 꼴로 다시 계산해 같은 값이 나오는지 본다.
prior_odds = prior_disease / (1 - prior_disease)
lr = sensitivity / false_positive_rate          # 가능도비
post_odds = prior_odds * lr
print(f"prior odds = {prior_odds:.6f},  LR = {lr:.2f},  post odds = {post_odds:.6f}")
print(f"posterior from odds = {post_odds / (1 + post_odds):.6f}")

# (2) 민감도를 1 로 올리면 가능도비가 10 이 된다. 정리 3의 상한이다.
best_odds = prior_odds * (1.0 / false_positive_rate)
print(f"sensitivity -> 1.0 : LR = {1.0/false_positive_rate:.2f},"
      f"  posterior = {best_odds / (1 + best_odds):.6f}")

P(disease | positive) = 0.0876
Despite a 95% sensitive test, only 8.8% of positives truly have the disease.
prior odds = 0.010101,  LR = 9.50,  post odds = 0.095960
posterior from odds = 0.087558
sensitivity -> 1.0 : LR = 10.00,  posterior = 0.091743

import numpy as np
import matplotlib.pyplot as plt

def bayes_update_visualization():
    """검사 성능은 그대로 두고 **유병률만** 바꿔 가며 사후확률을 그린다.

    "같은 검사, 같은 양성 결과인데 사후확률이 왜 달라지는가"에 답하는 그림이다.
    답은 분모에 있다. 유병률이 낮을수록 건강한 사람이 만들어 내는
    위양성이 분모를 지배해 사후확률을 끌어내린다.
    """
    prevalences = np.linspace(0.001, 0.5, 200)   # 0.1%부터 50%까지
    sensitivity = 0.95      # 검사 성능은 끝까지 고정
    specificity = 0.90

    posteriors = []
    for prev in prevalences:
        # 분모: 전확률의 법칙
        p_pos = sensitivity * prev + (1 - specificity) * (1 - prev)
        # 분자: 질병이면서 양성
        post = (sensitivity * prev) / p_pos
        posteriors.append(post)
    posteriors = np.array(posteriors)

    # 곡선이 50% 선을 넘는 첫 격자점과 (1)의 해석적 임계값을 견준다.
    i = int(np.argmax(posteriors >= 0.5))
    print(f"격자가 50%를 넘는 첫 유병률 = {prevalences[i]:.5f}"
          f"  (사후확률 {posteriors[i]:.4f})")
    print(f"해석적 임계 유병률         = "
          f"{(1 - specificity) / (sensitivity + 1 - specificity):.5f}")
    print(f"격자 간격                  = {prevalences[1] - prevalences[0]:.5f}")

    # (2) 원점 근처의 기울기. 어림은 LR x 유병률 이다.
    exact_1pct = sensitivity * 0.01 / (sensitivity * 0.01 + (1 - specificity) * 0.99)
    print(f"유병률 1%: 정확 {exact_1pct:.5f},  어림 9.5 x 0.01 = "
          f"{sensitivity / (1 - specificity) * 0.01:.5f}")

    fig, ax = plt.subplots(figsize=(12, 4))
    ax.plot(prevalences * 100, np.array(posteriors) * 100, lw=2)
    ax.set_xlabel('Prevalence (%)')
    ax.set_ylabel('P(Disease | Positive) (%)')
    ax.set_title("Bayes' Theorem: Posterior vs. Prevalence")
    # 50% 기준선. 이 선을 넘어야 "양성이면 아플 가능성이 더 크다"고 말할 수 있다.
    ax.axhline(y=50, color='r', linestyle='--', alpha=0.5, label='50% threshold')
    ax.legend()
    ax.spines[['top', 'right']].set_visible(False)
    plt.tight_layout()
    plt.show()

bayes_update_visualization()

격자가 50%를 넘는 첫 유병률 = 0.09629  (사후확률 0.5030)
해석적 임계 유병률         = 0.09524
격자 간격                  = 0.00251
유병률 1%: 정확 0.08756,  어림 9.5 x 0.01 = 0.09500

import numpy as np

p0, p1 = 0.5, 0.7                                # 공정한 동전 vs 치우친 동전

def log_bf(k, n):
    return k * np.log(p1 / p0) + (n - k) * np.log((1 - p1) / (1 - p0))

print("앞면 비율이 계속 0.7 로 관측될 때 (사전확률 0.5)")
print(f"{'n':>6}{'앞면':>6}{'로그 베이즈 인자':>18}{'사후확률':>14}")
for n in (10, 20, 50, 100, 200):
    k = int(round(0.7 * n))
    lbf = log_bf(k, n)
    print(f"{n:>6}{k:>6}{lbf:>18.4f}{1 / (1 + np.exp(-lbf)):>14.6f}")

print()
for target in (0.95, 0.99, 0.999):
    for n in range(1, 500):
        k = int(round(0.7 * n))
        if 1 / (1 + np.exp(-log_bf(k, n))) >= target:
            print(f"  사후확률 {target}에 이르려면 n = {n}")
            break

앞면 비율이 계속 0.7 로 관측될 때 (사전확률 0.5)
     n    앞면         로그 베이즈 인자          사후확률
    10     7            0.8228      0.694836
    20    14            1.6457      0.838303
    50    35            4.1141      0.983923
   100    70            8.2283      0.999733
   200   140           16.4566      1.000000

  사후확률 0.95에 이르려면 n = 34
  사후확률 0.99에 이르려면 n = 54
  사후확률 0.999에 이르려면 n = 81

import numpy as np

p0, p1 = 0.5, 0.7
log_bf = lambda k, n: k * np.log(p1 / p0) + (n - k) * np.log((1 - p1) / (1 - p0))

print("앞면 비율 0.7 이 관측될 때, 사전확률별 사후확률")
print(f"{'사전':>8}" + "".join(f"{'n=' + str(n):>13}" for n in (10, 30, 100, 300)))
for prior in (0.01, 0.10, 0.50, 0.90):
    row = f"{prior:>8.2f}"
    for n in (10, 30, 100, 300):
        k = int(round(0.7 * n))
        odds = prior / (1 - prior) * np.exp(log_bf(k, n))
        row += f"{odds / (1 + odds):>13.6f}"
    print(row)

앞면 비율 0.7 이 관측될 때, 사전확률별 사후확률
      사전         n=10         n=30        n=100        n=300
    0.01     0.022482     0.106535     0.974248     1.000000
    0.10     0.201911     0.567403     0.997603     1.000000
    0.50     0.694836     0.921903     0.999733     1.000000
    0.90     0.953472     0.990675     0.999970     1.000000

import numpy as np

prevalence, sens, spec = 0.001, 0.99, 0.95
ppv = lambda pr, se, sp: se * pr / (se * pr + (1 - sp) * (1 - pr))

p1 = ppv(prevalence, sens, spec)
p2 = ppv(p1, sens, spec)
p3 = ppv(p2, sens, spec)
print("검사를 반복하면 (오류가 서로 독립이라고 가정)")
print(f"  1차 양성 후 {p1:.6f}")
print(f"  2차도 양성  {p2:.6f}")
print(f"  3차도 양성  {p3:.6f}")

print("\n오류가 상관되어 있으면 2차의 조건부 특이도가 낮아진다")
for sp2 in (0.95, 0.70, 0.50, 0.20):
    print(f"  조건부 특이도 {sp2:.2f}: 사후확률 {ppv(p1, sens, sp2):.6f}")

print(f"\n{'유병률':>9}{'PPV':>10}{'NPV':>11}")
for pr in (0.0001, 0.001, 0.01, 0.05, 0.20):
    npv = spec * (1 - pr) / (spec * (1 - pr) + (1 - sens) * pr)
    print(f"{pr:>9.4f}{ppv(pr, sens, spec):>10.5f}{npv:>11.6f}")

검사를 반복하면 (오류가 서로 독립이라고 가정)
  1차 양성 후 0.019435
  2차도 양성  0.281832
  3차도 양성  0.885977

오류가 상관되어 있으면 2차의 조건부 특이도가 낮아진다
  조건부 특이도 0.95: 사후확률 0.281832
  조건부 특이도 0.70: 사후확률 0.061390
  조건부 특이도 0.50: 사후확률 0.037761
  조건부 특이도 0.20: 사후확률 0.023940

      유병률       PPV        NPV
   0.0001   0.00198   0.999999
   0.0010   0.01943   0.999989
   0.0100   0.16667   0.999894
   0.0500   0.51031   0.999446
   0.2000   0.83193   0.997375

import numpy as np
from scipy import stats

a, b = 2.0, 2.0                                  # 약한 사전: 평균 0.5
print(f"사전 Beta({a:.0f}, {b:.0f}): 평균 {a/(a+b):.4f}")
print(f"\n{'자료':>12}{'사후':>18}{'사후 평균':>12}{'95% 구간':>22}")
for k, n in [(7, 10), (70, 100), (700, 1000)]:
    A, B = a + k, b + n - k
    lo, hi = stats.beta.ppf([0.025, 0.975], A, B)
    print(f"{f'{k}/{n}':>12}{f'Beta({A:.0f}, {B:.0f})':>18}"
          f"{A/(A+B):>12.5f}   [{lo:.4f}, {hi:.4f}]")
print(f"\n표본비율은 세 경우 모두 0.700")

사전 Beta(2, 2): 평균 0.5000

          자료                사후       사후 평균                95% 구간
        7/10        Beta(9, 5)     0.64286   [0.3857, 0.8614]
      70/100      Beta(72, 32)     0.69231   [0.6008, 0.7769]
    700/1000    Beta(702, 302)     0.69920   [0.6705, 0.7272]

표본비율은 세 경우 모두 0.700