최소제곱 모의실험 (몬테카를로)¶

import numpy as np

def gen_X(n, k):
    """설계행렬. 첫 열의 1 이 절편에 대응한다."""
    return np.hstack([np.ones((n, 1)), np.random.randn(n, k - 1)])

def ols(y, X):
    """정규방정식의 해. 실제 계산에서는 역행렬 대신 solve 를 쓰는 편이 낫다."""
    return np.linalg.inv(X.T @ X) @ X.T @ y

def proj_P(X):
    """사영행렬 P. y 를 X 의 열공간 위로 떨어뜨린다. P @ y 가 곧 적합값이다."""
    return X @ np.linalg.inv(X.T @ X) @ X.T

def proj_M(X):
    """잔차생성행렬 M = I - P. M @ y 가 잔차이고, 열공간에 수직이다."""
    return np.eye(X.shape[0]) - proj_P(X)

def anova_decomposition(y, X, beta_hat):
    """TSS = ESS + RSS 로 갈라 본다.

    적합값과 잔차가 서로 수직이므로 피타고라스 정리가 그대로 성립한다.
    최소제곱의 기하가 이 한 줄에 들어 있다.
    """
    y_bar = y.mean()
    y_hat = X @ beta_hat
    TSS = float(np.sum((y - y_bar) ** 2))
    ESS = float(np.sum((y_hat - y_bar) ** 2))
    RSS = float(np.sum((y - y_hat) ** 2))
    return TSS, ESS, RSS

# P 와 M 의 성질, 그리고 분산분석 분해가 절편에 기대고 있음을 확인한다.
np.random.seed(0)
n, k = 50, 3
X = gen_X(n, k)
y = X @ np.array([[2.0], [3.0], [-1.0]]) + np.random.randn(n, 1)
P, M = proj_P(X), proj_M(X)
print(f"P 대칭  최대오차 {np.abs(P - P.T).max():.2e}     P 멱등  최대오차 {np.abs(P @ P - P).max():.2e}")
print(f"M 대칭  최대오차 {np.abs(M - M.T).max():.2e}     M 멱등  최대오차 {np.abs(M @ M - M).max():.2e}")
print(f"P @ M 의 최대 절대값 = {np.abs(P @ M).max():.2e}   (서로 직교하는 사영)")
print(f"tr(P) = {np.trace(P):.10f}  (k = {k}),   tr(M) = {np.trace(M):.10f}  (n-k = {n - k})")
print()
beta_hat = ols(y, X)
TSS, ESS, RSS = anova_decomposition(y, X, beta_hat)
print(f"절편 있음:  TSS = {TSS:.6f},  ESS + RSS = {ESS + RSS:.6f},"
      f"  차이 = {abs(TSS - (ESS + RSS)):.2e}")
e = y - X @ beta_hat
print(f"  잔차의 합 sum(e) = {e.sum():.2e}   (절편 열이 있으므로 0)")
print(f"  e 와 y-hat 의 내적 = {float((e.T @ (X @ beta_hat))[0, 0]):.2e}   (언제나 0)")
print()
X_noint = X[:, 1:]                     # 절편 열을 뺀다
beta_noint = ols(y, X_noint)
TSS2, ESS2, RSS2 = anova_decomposition(y, X_noint, beta_noint)
e2 = y - X_noint @ beta_noint
print(f"절편 없음:  TSS = {TSS2:.6f},  ESS + RSS = {ESS2 + RSS2:.6f},"
      f"  차이 = {TSS2 - (ESS2 + RSS2):+.6f}")
print(f"  잔차의 합 sum(e) = {e2.sum():.6f}   (0 이 아니다)")
print(f"  e 와 y-hat 의 내적 = {float((e2.T @ (X_noint @ beta_noint))[0, 0]):.2e}   (여전히 0)")
print(f"  교차항 -2*y-bar*sum(e) = {-2 * y.mean() * e2.sum():+.6f}   (차이와 같다)")

P 대칭  최대오차 2.78e-17     P 멱등  최대오차 8.33e-17
M 대칭  최대오차 2.78e-17     M 멱등  최대오차 5.55e-16
P @ M 의 최대 절대값 = 8.83e-17   (서로 직교하는 사영)
tr(P) = 3.0000000000  (k = 3),   tr(M) = 47.0000000000  (n-k = 47)

절편 있음:  TSS = 594.506762,  ESS + RSS = 594.506762,  차이 = 2.27e-13
  잔차의 합 sum(e) = 3.06e-14   (절편 열이 있으므로 0)
  e 와 y-hat 의 내적 = 1.49e-13   (언제나 0)

절편 없음:  TSS = 594.506762,  ESS + RSS = 1088.448126,  차이 = -493.941364
  잔차의 합 sum(e) = 112.048370   (0 이 아니다)
  e 와 y-hat 의 내적 = -4.29e-14   (여전히 0)
  교차항 -2*y-bar*sum(e) = -493.941364   (차이와 같다)

def monte_carlo(n=100, beta_true=[2, 3, -1], sigma=1.0, n_sim=5000):
    """같은 실험을 5000번 되풀이해 추정량의 분포를 본다.

    참 계수를 우리가 정해 두었으므로, 추정값들의 평균이 참값에 붙는지
    (불편성) 그리고 그 흩어짐이 얼마인지를 직접 확인할 수 있다.
    """
    k = len(beta_true)
    estimates = np.empty((n_sim, k))
    for i in range(n_sim):
        X = gen_X(n, k)
        beta = np.array(beta_true).reshape(-1, 1)
        u = np.random.randn(n, 1) * sigma
        y = X @ beta + u
        bhat = ols(y, X)
        estimates[i] = bhat.flatten()
    return estimates

beta_true = [2, 3, -1]
estimates = monte_carlo(n=200, beta_true=beta_true, sigma=2.0, n_sim=5000)
mc_mean = estimates.mean(axis=0)
mc_std = estimates.std(axis=0, ddof=1)

for j in range(len(beta_true)):
    print(f"beta_{j}: true={beta_true[j]}, "
          f"MC mean={mc_mean[j]:.4f}, MC std={mc_std[j]:.4f}")

beta_0: true=2, MC mean=2.0009, MC std=0.1414
beta_1: true=3, MC mean=3.0019, MC std=0.1425
beta_2: true=-1, MC mean=-0.9994, MC std=0.1442

# 이론값: 설명변수가 표준정규이므로 (X'X)^{-1} 의 대각원소의 기댓값이 알려져 있다.
n, k, sigma = 200, 3, 2.0
print(f"거친 근사   sigma/sqrt(n)     = {sigma / np.sqrt(n):.6f}")
print(f"기울기 둘   sigma/sqrt(n - 4) = {sigma / np.sqrt(n - 4):.6f}   "
      f"(E[(X'X)^-1_jj] = 1/(n-4))")
exact_00 = 1 / n + (k - 1) / (n * (n - 4))
print(f"절편        sigma*sqrt(1/n + (k-1)/(n(n-4))) = {sigma * np.sqrt(exact_00):.6f}")
print()
# 씨앗을 고정해 2000 번 되풀이한다 (쪽의 코드는 씨앗이 없어 재현되지 않는다)
np.random.seed(20251002)
reps = 2000
est = monte_carlo(n=n, beta_true=[2, 3, -1], sigma=sigma, n_sim=reps)
m, sd = est.mean(axis=0), est.std(axis=0, ddof=1)
theory = [sigma * np.sqrt(exact_00), sigma / np.sqrt(n - 4), sigma / np.sqrt(n - 4)]
truth = [2, 3, -1]
print(f"씨앗 고정, {reps} 회")
print(f"{'j':>3s}{'참값':>7s}{'MC 평균':>12s}{'MC 오차':>11s}{'z':>8s}"
      f"{'MC 표준편차':>14s}{'이론 SE':>11s}{'비':>9s}")
for j in range(k):
    mc_err = sd[j] / np.sqrt(reps)
    print(f"{j:3d}{truth[j]:7d}{m[j]:12.5f}{mc_err:11.5f}{(m[j] - truth[j]) / mc_err:8.2f}"
          f"{sd[j]:14.5f}{theory[j]:11.5f}{sd[j] / theory[j]:9.4f}")
print(f"  표준편차의 몬테카를로 상대오차 ~ 1/sqrt(2*{reps - 1}) = {1 / np.sqrt(2 * (reps - 1)):.4f}")
print()
# 씨앗이 없으면 같은 코드가 다른 수를 준다
a = monte_carlo(n=n, beta_true=[2, 3, -1], sigma=sigma, n_sim=200).mean(axis=0)
b = monte_carlo(n=n, beta_true=[2, 3, -1], sigma=sigma, n_sim=200).mean(axis=0)
print(f"씨앗 없이 200 회씩 두 번:  {np.round(a, 4)}  대  {np.round(b, 4)}")
print(f"  같은가? {np.allclose(a, b)}")

거친 근사   sigma/sqrt(n)     = 0.141421
기울기 둘   sigma/sqrt(n - 4) = 0.142857   (E[(X'X)^-1_jj] = 1/(n-4))
절편        sigma*sqrt(1/n + (k-1)/(n(n-4))) = 0.142141

씨앗 고정, 2000 회
  j     참값       MC 평균      MC 오차       z       MC 표준편차      이론 SE        비
  0      2     2.00224    0.00319    0.70       0.14249    0.14214   1.0025
  1      3     3.00504    0.00321    1.57       0.14338    0.14286   1.0036
  2     -1    -1.00354    0.00324   -1.09       0.14511    0.14286   1.0158
  표준편차의 몬테카를로 상대오차 ~ 1/sqrt(2*1999) = 0.0158

씨앗 없이 200 회씩 두 번:  [ 2.0101  2.9817 -0.9961]  대  [ 2.0039  2.9969 -1.0022]
  같은가? False

X = gen_X(50, 3)
P = proj_P(X)
M = proj_M(X)
print("P idempotent:", np.allclose(P @ P, P))
print("M idempotent:", np.allclose(M @ M, M))
print("P symmetric:", np.allclose(P, P.T))
print("M symmetric:", np.allclose(M, M.T))
print("tr(P):", np.trace(P))  # should be 3
print("tr(M):", np.trace(M))  # should be 47

P idempotent: True
M idempotent: True
P symmetric: True
M symmetric: True
tr(P): 3.000000000000001
tr(M): 46.99999999999999

from scipy import stats
coverage = np.zeros(3)
n, sigma = 200, 2.0
beta_true_arr = np.array(beta_true)
for i in range(5000):
    X = gen_X(n, 3)
    y = X @ beta_true_arr.reshape(-1, 1) + sigma * np.random.randn(n, 1)
    bhat = ols(y, X)
    e = y - X @ bhat
    s2 = np.sum(e ** 2) / (n - 3)
    se = np.sqrt(s2 * np.diag(np.linalg.inv(X.T @ X)))
    t_star = stats.t(n - 3).ppf(0.975)
    for j in range(3):
        if bhat[j, 0] - t_star * se[j] <= beta_true[j] <= bhat[j, 0] + t_star * se[j]:
            coverage[j] += 1
print("Coverage:", coverage / 5000)  # should be ~0.95

Coverage: [0.952  0.9482 0.951 ]