두 패러다임의 강점과 한계¶

"""같은 도구(선형모형)로 목표가 다른 두 일을 해 본다: 추론 대 예측."""

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n = 500
true_effect = 2.0

# === 고전적 A/B 검정 — 목표는 "효과가 있는가"를 판정하는 것 ===
# 처리를 무작위로 배정한다. 배정이 결과와 무관하므로 인과 해석이 가능하다.
group = rng.choice([0, 1], size=n)
outcome = 10 + true_effect * group + rng.normal(0, 5, n)

# 관심사는 계수 하나(효과의 크기)와 그것에 대한 불확실성(p값)이다.
# 예측 정확도는 아예 재지도 않는다. 잡음이 커서 개별 예측은 형편없다.
t_stat, p_val = stats.ttest_ind(outcome[group == 1], outcome[group == 0])
print("=== Classical A/B test ===")
print(f"Estimated effect: {outcome[group == 1].mean() - outcome[group == 0].mean():+.2f} "
      f"(true {true_effect})")
print(f"p = {p_val:.4f}")

# === 현대적 예측 — 목표는 "새 자료에서 얼마나 잘 맞히는가" ===
# 설명변수 3개짜리 회귀. 계수의 의미나 유의성에는 관심이 없다.
X = rng.standard_normal((n, 3))
y = 2 * X[:, 0] - X[:, 1] + 0.5 * X[:, 2] + rng.standard_normal(n)

# 자료를 훈련용과 시험용으로 반씩 나눈다.
# 이 분할이 예측 패러다임의 핵심이다. 모형이 **보지 않은** 자료로 평가해야
# 외운 것인지 배운 것인지 구별할 수 있다.
train, test = np.arange(n // 2), np.arange(n // 2, n)

# 훈련자료로만 계수를 구한다(최소제곱)
beta_train = np.linalg.lstsq(X[train], y[train], rcond=None)[0]

# 시험자료에서 평균제곱오차를 잰다. 이것이 성적표다.
mse_test = np.mean((y[test] - X[test] @ beta_train) ** 2)
print("\n=== Modern prediction ===")
print(f"Test MSE = {mse_test:.3f}")

# --- (1)(2) 의 이론값 ---
n1, n0 = int((group == 1).sum()), int((group == 0).sum())
print(f"\n군 크기 {n1} / {n0},  이론 SE = 5*sqrt(1/n1+1/n0) = "
      f"{5 * np.sqrt(1 / n1 + 1 / n0):.4f}")
p_dim, n_tr = X.shape[1], len(train)
print(f"이론 E[시험 MSE] = 1 + p/(n-p-1) = {1 + p_dim / (n_tr - p_dim - 1):.4f}")

=== Classical A/B test ===
Estimated effect: +2.35 (true 2.0)
p = 0.0000

=== Modern prediction ===
Test MSE = 1.017

군 크기 244 / 256,  이론 SE = 5*sqrt(1/n1+1/n0) = 0.4473
이론 E[시험 MSE] = 1 + p/(n-p-1) = 1.0122

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

INK, BLUE, ORANGE, GREEN, MUTED = "#37474F", "#1565C0", "#E65100", "#33691E", "#90A4AE"

rng = np.random.default_rng(0)
n = 60
x = np.sort(rng.uniform(-1, 1, n))
y = np.sin(3 * x) + rng.normal(0, 0.3, n)
xt = np.sort(rng.uniform(-1, 1, 2000))
yt = np.sin(3 * xt) + rng.normal(0, 0.3, 2000)

degrees = [1, 3, 5, 9, 15, 25]
train_mse, test_mse = [], []
print(f"{'다항식 차수':>11}{'훈련 MSE':>11}{'시험 MSE':>11}")
for d in degrees:
    c = np.polyfit(x, y, d)
    train_mse.append(np.mean((np.polyval(c, x) - y) ** 2))
    test_mse.append(np.mean((np.polyval(c, xt) - yt) ** 2))
    print(f"{d:>11}{train_mse[-1]:>11.4f}{test_mse[-1]:>11.4f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(degrees, train_mse, "o-", color=BLUE, label="훈련 MSE")
axes[0].plot(degrees, test_mse, "s-", color=ORANGE, label="시험 MSE")
axes[0].axhline(0.09, color=MUTED, ls=":", label="잡음 하한 $\\sigma^2 = 0.09$")
axes[0].set_xlabel("다항식 차수 (복잡도)")
axes[0].set_ylabel("평균제곱오차")
axes[0].set_title("훈련은 계속 내려가고 시험은 U자를 그린다", fontsize=10)
axes[0].legend(fontsize=8)

grid = np.linspace(-1, 1, 400)
for d, style, color in [(3, "-", BLUE), (25, "--", ORANGE)]:
    axes[1].plot(grid, np.polyval(np.polyfit(x, y, d), grid), style,
                 color=color, label=f"차수 {d}")
axes[1].scatter(x, y, s=18, color=INK, alpha=0.6, label="자료")
axes[1].plot(grid, np.sin(3 * grid), color=GREEN, lw=1.5, label="참 함수")
axes[1].set_ylim(-2, 2)
axes[1].set_xlabel("$x$")
axes[1].set_ylabel("$y$")
axes[1].set_title("차수 25 는 잡음까지 따라간다", fontsize=10)
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()

     다항식 차수     훈련 MSE     시험 MSE
          1     0.2318     0.2715
          3     0.0958     0.0979
          5     0.0837     0.0996
          9     0.0793     0.1025
         15     0.0744     0.1076
         25     0.0526     0.2551

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng(0)
K, n = 1000, 2000
true_effect = np.zeros(K)
true_effect[:50] = 0.15                 # 1000건 중 50건만 진짜 효과가 있다

pvals = np.empty(K)
for i in range(K):
    a = rng.normal(0, 1, n)
    b = rng.normal(true_effect[i], 1, n)
    pvals[i] = stats.ttest_ind(b, a)[1]

is_null = true_effect == 0
for label, sel in [
        ("보정 없음", pvals < 0.05),
        ("본페로니", multipletests(pvals, method="bonferroni")[0]),
        ("BH (FDR 0.05)", multipletests(pvals, method="fdr_bh")[0])]:
    n_sig = sel.sum()
    n_false = np.sum(sel & is_null)
    fdr = n_false / max(n_sig, 1)
    print(f"{label:>14}: 유의 {n_sig:>3}건   거짓 {n_false:>2}건   "
          f"실제 거짓발견율 {fdr:.3f}   참 효과 검출 {np.sum(sel & ~is_null):>2}/50")

         보정 없음: 유의 105건   거짓 55건   실제 거짓발견율 0.524   참 효과 검출 50/50
          본페로니: 유의  37건   거짓  0건   실제 거짓발견율 0.000   참 효과 검출 37/50
 BH (FDR 0.05): 유의  50건   거짓  3건   실제 거짓발견율 0.060   참 효과 검출 47/50

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold

rng = np.random.default_rng(0)
n, p, tau = 2000, 10, 1.0
rf = lambda: RandomForestRegressor(n_estimators=60, random_state=0, n_jobs=-1)

naive, plugin, dml = [], [], []
for _ in range(200):
    X = rng.normal(0, 1, (n, p))
    g = np.sin(X[:, 0]) + 0.5 * X[:, 1] ** 2 + X[:, 2]
    D = 0.8 * np.tanh(X[:, 0]) + 0.5 * X[:, 1] + rng.normal(0, 1, n)
    Y = tau * D + g + rng.normal(0, 1, n)

    naive.append(LinearRegression().fit(D.reshape(-1, 1), Y).coef_[0])

    # 순진한 플러그인: 같은 자료로 g 를 적합해 빼기
    gm = rf().fit(X, Y)
    plugin.append(LinearRegression().fit(
        D.reshape(-1, 1), Y - gm.predict(X)).coef_[0])

    # DML: 양쪽을 직교화하고 교차적합한다
    Yr, Dr = np.empty(n), np.empty(n)
    for tr, te in KFold(5, shuffle=True, random_state=0).split(X):
        Yr[te] = Y[te] - rf().fit(X[tr], Y[tr]).predict(X[te])
        Dr[te] = D[te] - rf().fit(X[tr], D[tr]).predict(X[te])
    dml.append((Dr @ Yr) / (Dr @ Dr))

for label, v in [("보정 없는 OLS", naive), ("순진한 플러그인", plugin),
                 ("DML (직교화 + 교차적합)", dml)]:
    v = np.array(v)
    print(f"{label:>22}: 평균 {v.mean():.4f}   편향 {v.mean() - tau:+.4f}"
          f"   표준편차 {v.std():.4f}")

             보정 없는 OLS: 평균 1.2173   편향 +0.2173   표준편차 0.0303
              순진한 플러그인: 평균 0.2514   편향 -0.7486   표준편차 0.0089
      DML (직교화 + 교차적합): 평균 0.9608   편향 -0.0392   표준편차 0.0269

import numpy as np
from sklearn.ensemble import RandomForestRegressor

rng = np.random.default_rng(1)
n = 20_000

X1 = rng.normal(0, 1, n)                 # 참 원인
U = rng.normal(0, 1, n)                  # 측정되지 않은 교란요인
X3 = U + rng.normal(0, 0.3, n)           # 교란요인의 대리
Y = 2.0 * X1 + 1.5 * U + rng.normal(0, 1, n)
X2 = Y + rng.normal(0, 0.5, n)           # Y 의 결과 — 인과효과는 정확히 0

F = np.column_stack([X1, X2, X3])
rf = RandomForestRegressor(n_estimators=200, random_state=0, n_jobs=-1).fit(F, Y)

names = ["X1  참 원인 (효과 2.0)", "X2  Y 의 결과 (효과 0)", "X3  교란 대리 (효과 0)"]
print("랜덤포레스트 변수 중요도")
for name, imp in sorted(zip(names, rf.feature_importances_), key=lambda z: -z[1]):
    print(f"  {name:<26} {imp:.4f}")
print(f"\n모형 결정계수 {rf.score(F, Y):.4f}   ← 예측은 거의 완벽하다")

랜덤포레스트 변수 중요도
  X2  Y 의 결과 (효과 0)          0.9763
  X1  참 원인 (효과 2.0)          0.0120
  X3  교란 대리 (효과 0)           0.0117

모형 결정계수 0.9956   ← 예측은 거의 완벽하다