조건부 독립¶

import numpy as np

def berkson_paradox_simulation(n_simulations=200_000):
    """버크슨의 역설: 독립인 두 사건이 조건을 걸면 종속이 된다.

    공통의 결과를 조건으로 삼으면, 원인 하나가 일어났다는 사실이
    나머지 원인의 확률을 낮춘다. 한쪽이 이미 설명해 버리기 때문이다.
    """
    np.random.seed(42)

    # 두 원인과 하나의 공통 결과.
    #   fire  = 화재 (드묾)
    #   toast = 토스트 태움 (흔함)
    #   alarm = 경보기 울림 (둘 중 하나만 일어나도 울린다)
    p_fire = 0.01
    p_toast = 0.10

    # 두 원인을 서로 **완전히 독립으로** 생성한다는 점이 중요하다.
    # 각자 따로 난수를 뽑으므로 설계상 아무 관계가 없다.
    fire = np.random.rand(n_simulations) < p_fire
    toast = np.random.rand(n_simulations) < p_toast
    alarm = fire | toast      # 공통 결과(충돌부, collider)

    # 조건 없이 보면 정말 독립이다.
    # fire[toast] 는 "토스트를 태운 시행들만" 골라 낸 것이고,
    # 그 안에서의 화재 비율이 전체 화재 비율과 같아야 한다.
    p_fire_given_toast = fire[toast].mean()
    print(f"P(fire) = {fire.mean():.4f}")
    print(f"P(fire | toast) = {p_fire_given_toast:.4f}")
    print(f"Unconditionally independent: {abs(fire.mean() - p_fire_given_toast) < 0.005}\n")

    # 그런데 공통 결과인 alarm 으로 조건을 걸면 둘이 얽힌다.
    # 경보가 울렸는데 토스트를 태우지 않았다면, 남은 설명은 화재뿐이다.
    # 한 원인을 배제하면 다른 원인의 확률이 치솟는 이 현상을
    # "설명해 없애기(explaining away)" 또는 버크슨의 역설이라 한다.
    p_fire_given_alarm = fire[alarm].mean()
    p_fire_given_alarm_no_toast = fire[alarm & ~toast].mean()
    print(f"P(fire | alarm) = {p_fire_given_alarm:.4f}")
    print(f"P(fire | alarm, no toast) = {p_fire_given_alarm_no_toast:.4f}")
    print(f"Conditionally dependent (explaining away): "
          f"{abs(p_fire_given_alarm - p_fire_given_alarm_no_toast) > 0.01}")

berkson_paradox_simulation()

P(fire) = 0.0098
P(fire | toast) = 0.0101
Unconditionally independent: True

P(fire | alarm) = 0.0899
P(fire | alarm, no toast) = 1.0000
Conditionally dependent (explaining away): True

import numpy as np

def mixture_coin_simulation(n_simulations=200_000):
    """혼합모형에서 조건부독립을 보인다.

    앞 절의 버크슨 역설과 정확히 반대 방향의 예다.
      버크슨: 독립이던 것이 조건을 걸면 종속이 된다 (공통 결과로 조건)
      여기  : 종속이던 것이 조건을 걸면 독립이 된다 (공통 원인으로 조건)
    """
    np.random.seed(42)

    # 숨은 공통 원인: 어느 동전을 골랐는가.
    #   동전 0은 앞면 확률 0.3, 동전 1은 0.7.
    # 관측자는 이 값을 볼 수 없고 던진 결과만 본다.
    coin = np.random.randint(0, 2, size=n_simulations)
    p_heads = np.where(coin == 0, 0.3, 0.7)

    # 같은 동전을 두 번 던진다. 두 번의 던짐은 동전이 정해지면 서로 무관하다.
    flip1 = np.random.rand(n_simulations) < p_heads
    flip2 = np.random.rand(n_simulations) < p_heads

    # 조건 없이 보면 종속이다.
    # 첫 번째가 앞면이면 "0.7짜리 동전일 가능성"이 커지고,
    # 그 정보가 두 번째 던짐의 예측을 바꾸기 때문이다.
    p_f2 = flip2.mean()
    p_f2_given_f1 = flip2[flip1].mean()
    print("=== Unconditional (marginal) ===")
    print(f"P(flip2=H) = {p_f2:.4f}")
    print(f"P(flip2=H | flip1=H) = {p_f2_given_f1:.4f}")
    print(f"Not independent: {abs(p_f2 - p_f2_given_f1) > 0.01}\n")

    # 그런데 동전이 무엇인지 알고 나면(= 공통 원인으로 조건을 걸면)
    # 첫 번째 결과가 더 알려 줄 것이 없어져 두 던짐이 독립이 된다.
    mask_c0 = coin == 0
    p_f2_c0 = flip2[mask_c0].mean()
    p_f2_given_f1_c0 = flip2[mask_c0 & flip1].mean()
    print("=== Conditional on coin 0 (P(H)=0.3) ===")
    print(f"P(flip2=H | coin=0) = {p_f2_c0:.4f}")
    print(f"P(flip2=H | flip1=H, coin=0) = {p_f2_given_f1_c0:.4f}")
    print(f"Conditionally independent: {abs(p_f2_c0 - p_f2_given_f1_c0) < 0.02}")

mixture_coin_simulation()

=== Unconditional (marginal) ===
P(flip2=H) = 0.5008
P(flip2=H | flip1=H) = 0.5809
Not independent: True

=== Conditional on coin 0 (P(H)=0.3) ===
P(flip2=H | coin=0) = 0.3007
P(flip2=H | flip1=H, coin=0) = 0.3004
Conditionally independent: True

import numpy as np

rng = np.random.default_rng(0)
N = 600_000

def partial_corr(a, b, c):
    """c 에 회귀한 잔차들의 상관 = c 를 조건으로 한 상관"""
    ra = a - np.polyval(np.polyfit(c, a, 1), c)
    rb = b - np.polyval(np.polyfit(c, b, 1), c)
    return np.corrcoef(ra, rb)[0, 1]

# 사슬 X → Z → Y
x1 = rng.normal(size=N); z1 = x1 + rng.normal(size=N); y1 = z1 + rng.normal(size=N)
# 포크 X ← Z → Y
z2 = rng.normal(size=N); x2 = z2 + rng.normal(size=N); y2 = z2 + rng.normal(size=N)
# 충돌부 X → Z ← Y
x3 = rng.normal(size=N); y3 = rng.normal(size=N); z3 = x3 + y3 + rng.normal(size=N)

print(f"{'구조':>18}{'주변 상관':>12}{'Z 조건부':>12}")
for name, (a, b, c) in [("사슬  X→Z→Y", (x1, y1, z1)),
                        ("포크  X←Z→Y", (x2, y2, z2)),
                        ("충돌부 X→Z←Y", (x3, y3, z3))]:
    print(f"{name:>18}{np.corrcoef(a, b)[0, 1]:>+12.4f}{partial_corr(a, b, c):>+12.4f}")

                구조       주변 상관       Z 조건부
         사슬  X→Z→Y     +0.5782     -0.0006
         포크  X←Z→Y     +0.5004     +0.0009
         충돌부 X→Z←Y     +0.0015     -0.4999

import numpy as np

rng = np.random.default_rng(2)
N = 600_000

def partial_corr(a, b, c):
    ra = a - np.polyval(np.polyfit(c, a, 1), c)
    rb = b - np.polyval(np.polyfit(c, b, 1), c)
    return np.corrcoef(ra, rb)[0, 1]

cases = {}
# (1) 둘 다 독립: X, Y, Z 모두 무관
x = rng.normal(size=N); y = rng.normal(size=N); z = rng.normal(size=N)
cases["주변 O, 조건부 O"] = (x, y, z)
# (2) 주변 독립, 조건부 종속: 충돌부
x = rng.normal(size=N); y = rng.normal(size=N); z = x + y + rng.normal(size=N)
cases["주변 O, 조건부 X"] = (x, y, z)
# (3) 주변 종속, 조건부 독립: 포크
z = rng.normal(size=N); x = z + rng.normal(size=N); y = z + rng.normal(size=N)
cases["주변 X, 조건부 O"] = (x, y, z)
# (4) 둘 다 종속: 포크 + 직접 경로
z = rng.normal(size=N); x = z + rng.normal(size=N); y = z + 0.8 * x + rng.normal(size=N)
cases["주변 X, 조건부 X"] = (x, y, z)

print(f"{'경우':>18}{'주변 상관':>12}{'조건부 상관':>14}")
for name, (a, b, c) in cases.items():
    print(f"{name:>18}{np.corrcoef(a, b)[0, 1]:>+12.4f}{partial_corr(a, b, c):>+14.4f}")

                경우       주변 상관        조건부 상관
       주변 O, 조건부 O     -0.0001       -0.0000
       주변 O, 조건부 X     +0.0000       -0.5003
       주변 X, 조건부 O     +0.5000       +0.0009
       주변 X, 조건부 X     +0.8328       +0.6254

import numpy as np
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(1)
n = 4000

print(f"{'클래스 내 상관':>14}{'나이브 베이즈':>15}{'LDA (상관 반영)':>18}{'확률 보정 오차':>16}")
for rho in (0.0, 0.5, 0.9, 0.99):
    y = rng.integers(0, 2, n)
    L = np.linalg.cholesky([[1, rho], [rho, 1]])
    X = rng.normal(size=(n, 2)) @ L.T + np.column_stack([y * 1.2, y * 1.2])

    acc_nb = cross_val_score(GaussianNB(), X, y, cv=5).mean()
    acc_lda = cross_val_score(LinearDiscriminantAnalysis(), X, y, cv=5).mean()

    proba = GaussianNB().fit(X, y).predict_proba(X)[:, 1]
    err = cnt = 0.0
    for lo in np.arange(0, 1, 0.2):
        m = (proba >= lo) & (proba < lo + 0.2)
        if m.sum() > 30:
            err += abs(proba[m].mean() - y[m].mean()) * m.sum()
            cnt += m.sum()
    print(f"{rho:>14.2f}{acc_nb:>15.4f}{acc_lda:>18.4f}{err / max(cnt, 1):>16.4f}")

      클래스 내 상관        나이브 베이즈       LDA (상관 반영)        확률 보정 오차
          0.00         0.8105            0.8095          0.0068
          0.50         0.7525            0.7550          0.0589
          0.90         0.7335            0.7343          0.0879
          0.99         0.7233            0.7238          0.1035

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
N = 400_000
n = 6

# 동전의 앞면 확률 theta 를 Beta(2,2) 에서 뽑고, 그 동전을 n 번 던진다
theta = rng.beta(2, 2, N)
X = (rng.random((N, n)) < theta[:, None]).astype(int)

print("교환가능성 확인: 순서만 다른 수열의 확률")
for pattern in [(1, 1, 0, 0, 0, 0), (0, 0, 1, 1, 0, 0), (0, 0, 0, 0, 1, 1)]:
    print(f"  {pattern}: {np.mean((X == np.array(pattern)).all(axis=1)):.6f}")

print(f"\n주변 독립인가?")
print(f"  P(X1=1)        = {X[:, 0].mean():.6f}")
print(f"  P(X1=1 | X2=1) = {X[X[:, 1] == 1, 0].mean():.6f}   ← 다르다")

print(f"\ntheta 를 조건으로 하면?")
for lo in (0.2, 0.5, 0.8):
    m = np.abs(theta - lo) < 0.02
    p1 = X[m, 0].mean()
    p1_given = X[m & (X[:, 1] == 1), 0].mean()
    print(f"  theta≈{lo}: P(X1=1)={p1:.4f}, P(X1=1|X2=1)={p1_given:.4f}   ← 같다")

교환가능성 확인: 순서만 다른 수열의 확률
  (1, 1, 0, 0, 0, 0): 0.011795
  (0, 0, 1, 1, 0, 0): 0.012190
  (0, 0, 0, 0, 1, 1): 0.012105

주변 독립인가?
  P(X1=1)        = 0.498875
  P(X1=1 | X2=1) = 0.598010   ← 다르다

theta 를 조건으로 하면?
  theta≈0.2: P(X1=1)=0.2038, P(X1=1|X2=1)=0.2044   ← 같다
  theta≈0.5: P(X1=1)=0.4969, P(X1=1|X2=1)=0.4962   ← 같다
  theta≈0.8: P(X1=1)=0.7959, P(X1=1|X2=1)=0.7984   ← 같다