조건부확률¶

import numpy as np

def conditional_probability_simulation(n_simulations=100_000):
    """주사위 두 개로 조건부확률을 모의실험한다.

    조건부확률을 계산하는 대신 **직접 세어 본다.**
    P(A|B) = (A와 B가 모두 일어난 횟수) / (B가 일어난 횟수)
    라는 정의를 그대로 코드로 옮긴 것이다.
    """
    np.random.seed(42)

    # randint(1, 7)은 1 이상 7 미만, 즉 1~6을 뽑는다
    die1 = np.random.randint(1, 7, size=n_simulations)
    die2 = np.random.randint(1, 7, size=n_simulations)
    total = die1 + die2

    # 여기가 조건부확률의 핵심이다. 표본공간을 "첫 주사위가 3인 시행"으로 좁힌다.
    mask_B = die1 == 3                            # 조건 B가 성립한 시행들
    mask_A_and_B = (die1 == 3) & (total == 8)     # 그중 A도 성립한 시행들

    # 전체 10만이 아니라 **좁혀진 표본공간의 크기**로 나눈다.
    # 이 한 줄이 P(A|B) = P(A∩B)/P(B) 를 실행한 것이다.
    p_conditional = mask_A_and_B.sum() / mask_B.sum()
    print(f"Simulated P(sum=8 | die1=3) = {p_conditional:.4f}")
    print(f"Theoretical P(sum=8 | die1=3) = {1/6:.4f}")

    # 분모가 몇인지, 그래서 오차가 얼마인지 함께 본다.
    n_B, n_AB = int(mask_B.sum()), int(mask_A_and_B.sum())
    se = np.sqrt((1/6) * (5/6) / n_B)
    print(f"n_B = {n_B:,} (기댓값 {n_simulations/6:,.0f}),  n_AB = {n_AB:,}")
    print(f"se = {se:.5f},  차이 = {p_conditional - 1/6:+.5f}"
          f"  ({(p_conditional - 1/6) / se:+.2f} se)")

conditional_probability_simulation()

Simulated P(sum=8 | die1=3) = 0.1705
Theoretical P(sum=8 | die1=3) = 0.1667
n_B = 16,390 (기댓값 16,667),  n_AB = 2,795
se = 0.00291,  차이 = +0.00386  (+1.33 se)

import numpy as np

def medical_test_simulation(n_people=1_000_000):
    """전확률의 법칙을 100만 명 모의실험으로 확인한다."""
    np.random.seed(42)

    prevalence = 0.01            # 유병률 P(질병) = 1%
    sensitivity = 0.95           # 민감도 P(양성|질병) = 95%
    false_positive_rate = 0.10   # 위양성률 P(양성|건강) = 10%

    # 1단계: 100만 명 각자가 질병을 가졌는지 정한다(확률 1%)
    has_disease = np.random.rand(n_people) < prevalence

    # 2단계: 검사 결과를 정한다. 여기서 확률이 사람에 따라 달라진다.
    # np.where(조건, 참일 때, 거짓일 때) 로 두 경우를 한 번에 처리한다.
    #   질병이 있으면 95% 확률로 양성
    #   건강하면    10% 확률로 양성
    test_positive = np.where(
        has_disease,
        np.random.rand(n_people) < sensitivity,
        np.random.rand(n_people) < false_positive_rate
    )

    # 전체 양성률. 전확률의 법칙이 예측하는 값과 맞는지 아래에서 비교한다.
    #   P(양성) = P(양성|질병)P(질병) + P(양성|건강)P(건강)
    #           = 0.95*0.01 + 0.10*0.99 = 0.1085
    p_positive = test_positive.mean()
    print(f"Simulated P(positive) = {p_positive:.4f}")
    print(f"Theoretical P(positive) = {0.1085:.4f}")

    # 합만이 아니라 두 경로를 따로 재어 본다.
    tp = (has_disease & test_positive).mean()        # 질병이면서 양성
    fp = (~has_disease & test_positive).mean()       # 건강하면서 양성
    se = np.sqrt(0.1085 * (1 - 0.1085) / n_people)
    print(f"  참양성   {tp:.5f} (이론 0.00950)")
    print(f"  거짓양성 {fp:.5f} (이론 0.09900)")
    print(f"  합       {tp + fp:.5f} (이론 0.10850),  se = {se:.5f},"
          f"  차이 = {p_positive - 0.1085:+.5f}")

medical_test_simulation()

Simulated P(positive) = 0.1085
Theoretical P(positive) = 0.1085
  참양성   0.00950 (이론 0.00950)
  거짓양성 0.09903 (이론 0.09900)
  합       0.10853 (이론 0.10850),  se = 0.00031,  차이 = +0.00004

import numpy as np

rng = np.random.default_rng(0)
N = 2_000_000
kids = rng.integers(0, 2, (N, 2))            # 0=딸, 1=아들
day = rng.integers(0, 7, (N, 2))             # 0~6, 1=화요일
both = kids.sum(axis=1) == 2

cond_a = kids.sum(axis=1) >= 1
cond_b = kids[:, 0] == 1
cond_c = ((kids == 1) & (day == 1)).any(axis=1)

print(f"(a) 적어도 하나 아들:        {both[cond_a].mean():.6f}   이론 1/3 = {1/3:.6f}")
print(f"(b) 첫째가 아들:             {both[cond_b].mean():.6f}   이론 1/2 = {0.5:.6f}")
print(f"(c) 화요일에 태어난 아들 있음: {both[cond_c].mean():.6f}   이론 13/27 = {13/27:.6f}")

(a) 적어도 하나 아들:        0.333428   이론 1/3 = 0.333333
(b) 첫째가 아들:             0.500179   이론 1/2 = 0.500000
(c) 화요일에 태어난 아들 있음: 0.481286   이론 13/27 = 0.481481

import numpy as np
from math import comb

# 카드 13장 하트 중 연속으로 뽑기
print("52장에서 비복원으로 k장을 뽑아 모두 하트일 확률")
print(f"{'k':>4}{'연쇄 법칙':>14}{'조합으로':>14}")
for k in (1, 2, 3, 5, 13):
    chain = np.prod([(13 - i) / (52 - i) for i in range(k)])
    direct = comb(13, k) / comb(52, k)
    print(f"{k:>4}{chain:>14.8f}{direct:>14.8f}")

# 마르코프 가정이 있으면 항이 줄어든다
print("\n연쇄 법칙의 항 개수 (조건에 들어가는 변수 수)")
print(f"{'n':>5}{'일반':>12}{'마르코프(1차)':>16}")
for n in (3, 5, 10, 50):
    print(f"{n:>5}{sum(range(n)):>12}{n - 1:>16}")

52장에서 비복원으로 k장을 뽑아 모두 하트일 확률
   k         연쇄 법칙          조합으로
   1    0.25000000    0.25000000
   2    0.05882353    0.05882353
   3    0.01294118    0.01294118
   5    0.00049520    0.00049520
  13    0.00000000    0.00000000

연쇄 법칙의 항 개수 (조건에 들어가는 변수 수)
    n          일반        마르코프(1차)
    3           3               2
    5          10               4
   10          45               9
   50        1225              49

import numpy as np

rng = np.random.default_rng(2)
N = 4_000_000
x = rng.normal(0, 1, N)
y = 0.8 * x + rng.normal(0, 0.6, N)

print(f"P(X = 0.5) 를 추정하려 하면")
for eps in (0.1, 0.01, 0.001, 0.0001):
    m = np.abs(x - 0.5) < eps
    print(f"  |X - 0.5| < {eps:<7}: 해당 표본 {m.sum():>8}개 "
          f"({m.mean():.6f})   E[Y|·] = {y[m].mean() if m.sum() else float('nan'):.4f}")
print(f"\n이론값 E[Y | X=0.5] = 0.8 x 0.5 = {0.8 * 0.5:.4f}")

P(X = 0.5) 를 추정하려 하면
  |X - 0.5| < 0.1    : 해당 표본   281692개 (0.070423)   E[Y|·] = 0.3988
  |X - 0.5| < 0.01   : 해당 표본    28047개 (0.007012)   E[Y|·] = 0.4044
  |X - 0.5| < 0.001  : 해당 표본     2814개 (0.000704)   E[Y|·] = 0.4032
  |X - 0.5| < 0.0001 : 해당 표본      256개 (0.000064)   E[Y|·] = 0.4033

이론값 E[Y | X=0.5] = 0.8 x 0.5 = 0.4000

import numpy as np

def theory(p, i, N):
    if abs(p - 0.5) < 1e-12:
        return i / N
    q = (1 - p) / p
    return (1 - q ** i) / (1 - q ** N)

def simulate(p, i, N, reps=40_000, seed=1):
    rng = np.random.default_rng(seed)
    wins = 0
    for _ in range(reps):
        x = i
        while 0 < x < N:
            x += 1 if rng.random() < p else -1
        wins += (x == N)
    return wins / reps

print(f"{'p':>6}{'시작':>6}{'목표':>6}{'모의':>10}{'이론':>10}")
for p, i, N in [(0.5, 5, 10), (0.45, 5, 10), (0.6, 5, 10)]:
    print(f"{p:>6}{i:>6}{N:>6}{simulate(p, i, N):>10.4f}{theory(p, i, N):>10.4f}")

print("\n공정한 게임에서 목표를 키우면")
for N in (10, 50, 200, 1000):
    print(f"  1원으로 시작해 {N}원 목표: {theory(0.5, 1, N):.6f}")

print("\n약간 불리한 게임(p=0.49)에서 자본을 늘리면")
for i, N in [(10, 20), (50, 100), (200, 400), (500, 1000)]:
    print(f"  {i}원으로 시작해 {N}원 목표: {theory(0.49, i, N):.3e}")

     p    시작    목표        모의        이론
   0.5     5    10    0.5017    0.5000
  0.45     5    10    0.2688    0.2683
   0.6     5    10    0.8845    0.8836

공정한 게임에서 목표를 키우면
  1원으로 시작해 10원 목표: 0.100000
  1원으로 시작해 50원 목표: 0.020000
  1원으로 시작해 200원 목표: 0.005000
  1원으로 시작해 1000원 목표: 0.001000

약간 불리한 게임(p=0.49)에서 자본을 늘리면
  10원으로 시작해 20원 목표: 4.013e-01
  50원으로 시작해 100원 목표: 1.192e-01
  200원으로 시작해 400원 목표: 3.350e-04
  500원으로 시작해 1000원 목표: 2.056e-09