1 이 절의 위치
§10.1에서는 MLE의 점근 효율성을 다뤘다 — 모형이 맞을 때 MLE가 가장 효율적이다. 그런데 모형이 틀렸다면 어떻게 되는가?
§10.2 강건성(robustness)은 이 질문에 답한다. 모형에서 약간 벗어났을 때도 합리적인 성능을 유지하는 추정량을 강건 추정량(robust estimator)이라 한다 (Casella & Berger, 2002, Ch.10).
2 강건성의 정의와 세 가지 기준
Huber (1981)는 훌륭한 통계적 절차가 갖춰야 할 세 가지 조건을 제시했다:
- 최적성: 가정된 모형에서 합리적으로 좋은(최적 또는 근사 최적) 효율을 가져야 한다.
- 안정성: 모형 가정에서 작은 이탈이 있어도 성능이 조금만 저하되어야 한다.
- 재난 방지: 모형에서 다소 큰 이탈이 있어도 파국적 결과를 초래하지 않아야 한다.
강건성과 최적성은 trade-off 관계이다. 강건한 추정량은 모형이 정확할 때 최적 추정량보다 열등하지만, 모형이 틀렸을 때 훨씬 더 안정적으로 작동한다. 어느 쪽이 더 중요한지는 문제마다 다르다.
3 표본 평균의 강건성 분석
3.1 기준 1: 가정된 모형에서의 최적성
\(X_1, \ldots, X_n \stackrel{\text{iid}}{\sim} N(\mu, \sigma^2)\) 이면 \(\bar X\) 는 크래머-라오 하한을 달성하는 UMVUE이다. 기준 1은 완벽하게 충족된다.
3.2 기준 2: ε-오염 모형 (작은 이탈)
모형 이탈을 정량화하는 표준적 방법은 ε-오염 모형(ε-contamination model)이다. 각 관측값이 주 분포에서 \((1-\delta)\) 확률로, 오염 분포 \(f(x)\) 에서 \(\delta\) 확률로 생성된다:
\[ X_i \sim \begin{cases} N(\mu, \sigma^2) & \text{확률 } 1-\delta \\ f(x) & \text{확률 } \delta \end{cases} \]
오염 분포 \(f\) 가 평균 \(\theta\), 분산 \(\tau^2\) 를 갖는다면, 혼합 분포에서 \(\bar X\) 의 분산은
\[ \text{Var}(\bar X) = (1-\delta)\frac{\sigma^2}{n} + \delta\frac{\tau^2}{n} + \frac{\delta(1-\delta)(\theta-\mu)^2}{n}. \]
\(\theta \approx \mu\), \(\tau \approx \sigma\) 이면 \(\bar X\) 는 여전히 좋은 추정량이다. 그러나 \(f\) 가 코시 분포(Cauchy distribution)이면 \(\text{Var}(X_i) = \infty\) 이므로 \(\text{Var}(\bar X) = \infty\) 가 된다 — 기준 2가 실패한다.
직관: 코시 분포는 꼬리가 매우 두꺼워 이상치가 극단적으로 크게 나올 수 있다. 단 한 번의 코시 표본이 섞여도 \(\bar X\) 의 분산이 무한대로 폭발한다.
3.3 기준 3: 붕괴점 (Breakdown Value)
“가장 큰 관측값 \(X_{(n)}\)이 무한대로 커지면 \(\bar X\) 는 어떻게 되는가?” — 마찬가지로 무한대가 된다. 이 현상을 정식으로 정의한 것이 붕괴점(breakdown value)이다.
순서통계량 \(X_{(1)} < \cdots < X_{(n)}\) 에 기반한 통계량 \(T_n\) 의 붕괴점 \(b\) \((0 \leq b \leq 1)\) 는:
\[ \lim_{X_{(\lfloor(1-b)n\rfloor)} \to \infty} T_n < \infty \quad \text{이고} \quad \lim_{X_{(\lfloor(1-(b+\varepsilon))n\rfloor)} \to \infty} T_n = \infty \]
를 만족하는 값이다.
즉, 표본의 \((1-b)\) 분위수가 무한대로 커져도 \(T_n\) 이 유한하게 유지되는 최대 비율이다.
직관: 붕괴점이 높을수록 더 많은 관측값이 망가져도 추정량이 살아남는다.
| 추정량 | 붕괴점 | 의미 |
|---|---|---|
| 표본 평균 \(\bar X\) | 0 | 단 1개의 관측값만 무한대로 가도 \(\bar X \to \infty\) |
| 표본 중앙값 \(M_n\) | 50% | 표본의 절반이 극단값이 되어도 중앙값은 움직이지 않음 |
| 절사평균 (α-trimmed mean) | α | 양쪽 α% 를 제거하고 평균을 냄 |
표본 평균의 붕괴점이 0이라는 것은 기준 3의 완전한 실패이다. 단 한 개의 “오염된” 관측값이 \(\bar X\) 를 원하는 방향으로 얼마든지 끌어당길 수 있다.
4 표본 중앙값의 점근 정규성
중앙값은 붕괴점이 50%이므로 강건성이 뛰어나다. 그렇다면 정규 모형에서 얼마나 효율성을 잃는가? 이를 계산하려면 중앙값의 점근 분포가 필요하다.
모집단 pdf \(f\), cdf \(F\) (미분가능), \(F(\mu) = 1/2\) (즉 \(\mu\) 는 모집단 중앙값)인 경우, 표본 중앙값 \(M_n\) 은 점근적으로 다음을 따른다:
\[ \sqrt{n}(M_n - \mu) \xrightarrow{d} N\!\left(0,\; \frac{1}{4f(\mu)^2}\right). \]
증명 스케치:
이항 근사 논증을 이용한다. \(Y_i = \mathbf{1}[X_i \leq \mu + a/\sqrt{n}]\) 로 정의하면 \(Y_i \sim \text{Bernoulli}(p_n)\), \(p_n = F(\mu + a/\sqrt{n}) \to 1/2\).
\(P(\sqrt{n}(M_n - \mu) \leq a) = P\!\left(\sum_i Y_i \geq \tfrac{n+1}{2}\right)\)
을 표준화한 후 CLT를 적용하면:
\[ \frac{(n+1)/2 - np_n}{\sqrt{np_n(1-p_n)}} \to -2af(\mu) \]
임을 보일 수 있다. 따라서
\[ P(\sqrt{n}(M_n - \mu) \leq a) \to P(Z \geq -2af(\mu)) = \Phi(2af(\mu)), \]
이는 \(N(0, 1/[4f(\mu)^2])\) 의 CDF이다.
해석: 중앙값의 점근 분산은 \(f(\mu)^{-2}/4\) 이다. 모집단 pdf 의 중앙값 높이 \(f(\mu)\) 가 클수록 중앙값 추정이 더 정확해진다 — 그 지점에서 분포가 “뾰족”하여 중앙값 위치가 잘 정해지기 때문이다.
5 평균 vs 중앙값: ARE 비교
표본 평균의 점근 분산: \(\sigma^2/n\), 즉 \(\sqrt{n}(\bar X - \mu) \to N(0, \sigma^2)\).
따라서 ARE(중앙값, 평균)는:
\[ \text{ARE}(M_n, \bar X) = \frac{\sigma^2}{1/[4f(\mu)^2]} = 4\sigma^2 f(\mu)^2. \]
세 가지 대표 분포에 대해 계산하면:
| 분포 | \(f(\mu)\) | \(\sigma^2\) | ARE(중앙값, 평균) |
|---|---|---|---|
| 정규 \(N(0,1)\) | \(\frac{1}{\sqrt{2\pi}}\) | 1 | \(\frac{2}{\pi} \approx 0.637\) |
| 로지스틱 | \(\frac{1}{4}\) | \(\frac{\pi^2}{3}\) | \(\frac{\pi^2}{12} \approx 0.822\) |
| 이중지수 (라플라스) | \(\frac{1}{2}\) | 2 | 2 |
해석:
- 정규 분포: 중앙값은 평균보다 약 36% 비효율적이다. 100개 표본으로 얻는 중앙값의 정밀도를 평균으로 얻으려면 약 64개면 충분하다.
- 로지스틱 분포: 꼬리가 정규보다 약간 두꺼워 중앙값의 효율성 손실이 작아진다.
- 이중지수 분포: 중앙값이 평균보다 2배 효율적이다. 이중지수 분포에서 중앙값이 MLE이기 때문이다.
일반 원리: 분포의 꼬리가 두꺼울수록 평균 대비 중앙값의 효율성이 높아진다.
6 M-추정량: 강건성과 효율성의 타협
6.1 통합 프레임워크
평균과 중앙값, MLE는 모두 손실 함수(loss function) \(\rho\) 의 최소화로 표현된다:
\[ \hat\theta = \arg\min_a \sum_{i=1}^n \rho(x_i - a). \]
| 추정량 | \(\rho(x)\) | 특징 |
|---|---|---|
| 표본 평균 | \(\frac{1}{2}x^2\) | 이상치에 민감 (제곱 패널티) |
| 표본 중앙값 | \(|x|\) | 이상치에 강건 (절댓값 패널티) |
| MLE | \(-\log f(x|\theta)\) | 모형이 맞으면 최적 |
제곱 손실은 큰 잔차에 과도한 가중치를 부여해 이상치에 취약하다. 절댓값 손실은 이상치를 무시하지만 미분 불가능하고 정규 모형에서 비효율적이다.
Huber (1964)는 이 둘의 타협점을 제안했다.
6.2 Huber 추정량
\[ \rho(x) = \begin{cases} \frac{1}{2}x^2 & \text{if } |x| \leq k \\ k|x| - \frac{1}{2}k^2 & \text{if } |x| > k \end{cases} \]
- \(|x| \leq k\) (작은 잔차): 제곱 손실처럼 작동 → 효율성 확보
- \(|x| > k\) (큰 잔차): 절댓값 손실처럼 작동 → 이상치 내성
- 튜닝 상수 \(k\): 작을수록 중앙값에 가깝고 (강건), 클수록 평균에 가깝다 (효율적)
최솟값 조건 \(\sum_i \psi(x_i - \theta) = 0\) 에서 \(\psi = \rho'\):
\[ \psi(x) = \begin{cases} x & \text{if } |x| \leq k \\ k \cdot \text{sign}(x) & \text{if } |x| > k \end{cases} \]
이 \(\psi\) 는 잔차가 \(k\) 를 초과하면 더 이상 패널티를 키우지 않고 \(\pm k\) 로 제한(clamp)한다. 따라서 단일 이상치가 미치는 영향이 \(|k|\) 이하로 bounded된다.
예시 (Ex 10.2.5): 데이터 = 8개 표준 정규 관측값 + 3개 이상치 \(\{3, 6, 9\}\)
\[ \mathbf{x} = \{-1.28, -0.96, -0.46, -0.44, -0.26, -0.21, -0.063, 0.39, 3, 6, 9\} \]
| \(k\) | 0 (중앙값) | 1 | 3 | 5 | 10 | \(\infty\) (평균) |
|---|---|---|---|---|---|---|
| 추정값 | -0.21 | 0.03 | 0.29 | 0.52 | 0.97 | 1.33 |
\(k\) 가 커질수록 이상치의 영향이 커져 평균(1.33)에 가까워진다. 이상치 없는 정규 관측값의 진짜 중심은 약 0 근처이므로, 작은 \(k\) 값이 더 정확하다.
7 M-추정량의 점근 정규성
일반 M-추정량 \(\hat\theta_M\) (손실 함수 \(\rho\), 기울기 함수 \(\psi = \rho'\)) 에 대해 Taylor 전개 논증(§10.1.2와 동일한 방식)을 적용하면:
\(E_{\theta_0}[\psi(X - \theta_0)] = 0\) 이 성립할 때,
\[ \sqrt{n}(\hat\theta_M - \theta_0) \xrightarrow{d} N\!\left(0,\; \frac{E_{\theta_0}[\psi(X-\theta_0)^2]}{[E_{\theta_0}\psi'(X-\theta_0)]^2}\right). \]
특수 케이스 확인:
- \(\psi(x) = x\) (평균): \(E[\psi^2] = \sigma^2\), \(E[\psi'] = 1\) → 분산 = \(\sigma^2\) ✓
- \(\psi(x) = \text{sign}(x)\) (중앙값): \(E[\psi^2] = 1\), \(E[\psi'] = 2f(\mu)\) → 분산 = \(1/(4f(\mu)^2)\) ✓
- \(\psi(x) = l'(\theta_0|x)\) (MLE): \(E[\psi^2] = I(\theta_0)\), \(E[\psi'] = -I(\theta_0)\) → 분산 = \(1/I(\theta_0)\) ✓
8 M-추정량 vs MLE: 효율성 상한
M-추정량의 ARE 대비 MLE는 항상 1 이하이다.
\[ \text{ARE}(\hat\theta_M, \hat\theta_{\text{MLE}}) = \frac{[E_\theta \psi(X-\theta_0)l'(\theta|X)]^2} {E_\theta\psi(X-\theta)^2 \cdot E_\theta[l'(\theta|X)]^2} \leq 1. \]
코시-슈바르츠 부등식 \([E(UV)]^2 \leq E(U^2)E(V^2)\) 에 의해 등호는 \(\psi \propto l'\) 일 때만 성립한다 — 즉, M-추정량의 기울기 함수가 점수 함수에 비례할 때.
Huber 추정량의 ARE (\(k = 1.5\)):
| 분포 | vs 평균 | vs 중앙값 |
|---|---|---|
| 정규 | 0.96 | 1.51 |
| 로지스틱 | 1.08 | 1.31 |
| 이중지수 | 1.37 | 0.68 |
Huber 추정량 (\(k=1.5\))은: - 정규 분포에서 평균과 거의 동등한 효율 (4% 손실만) - 로지스틱과 이중지수에서는 평균보다 더 효율적 - 이중지수에서 중앙값보다 덜 효율적 (중앙값이 MLE이므로 예상 가능)
이것이 Huber 추정량의 강점이다 — 어떤 분포에서도 합리적인 성능을 보장한다.
9 강건성 지표 요약
| 추정량 | 붕괴점 | 정규에서의 ARE | 코시 모형 |
|---|---|---|---|
| 표본 평균 | 0% | 1.00 (최적) | \(\text{Var} = \infty\) |
| 표본 중앙값 | 50% | 0.637 | 유한 |
| Huber (\(k=1.5\)) | ~7.5%* | 0.96 | 유한 |
| α-절사평균 | α | \(1 - 2\alpha\) 근방 | 유한 |
*Huber 붕괴점은 \(k\) 에 따라 다름. \(k \to 0\) 이면 50%, \(k \to \infty\) 이면 0%에 가까워짐.
10 응용 분야
| 분야 | 활용 | 이유 |
|---|---|---|
| 금융 데이터 분석 | 수익률 위치 추정 | 시장 충격으로 극단 수익률 빈번 |
| 의료 임상시험 | 혈압·혈당 효과 추정 | 오측값, 프로토콜 위반 관측값 |
| 이미지 처리 | 노이즈 제거 | 소금-후추 노이즈는 충격 잡음 |
| 회귀 분석 | 이상치 영향 최소화 | Huber 회귀, RANSAC |
| 기후 과학 | 기온 이상 탐지 | 센서 오류, 극단 기상 이벤트 |
11 코드 예시
11.1 Step 1: 순수 Python 구현 (원리 이해)
import numpy as np
# ε-오염 모형: 표본 평균 분산 폭발 시뮬레이션
def simulate_contamination(n=100, delta=0.05, n_sim=2000, seed=42):
"""
ε-오염 모형에서 표본 평균 vs 중앙값의 분산 비교.
Parameters
----------
n : 표본 크기
delta : 오염 비율 (Cauchy 관측값 비율)
n_sim : 시뮬레이션 반복 수
seed : 난수 시드
"""
rng = np.random.default_rng(seed)
means, medians = [], []
for _ in range(n_sim):
# 주 분포: N(0,1)
x_normal = rng.normal(0, 1, size=n)
# 오염: Cauchy(0,1) — 분산 무한대
n_contaminated = int(delta * n)
x_cauchy = rng.standard_cauchy(size=n_contaminated)
# 오염 표본: (1-δ) N(0,1) + δ Cauchy
idx = rng.choice(n, size=n_contaminated, replace=False)
x_contaminated = x_normal.copy()
x_contaminated[idx] = x_cauchy
means.append(np.mean(x_contaminated))
medians.append(np.median(x_contaminated))
# 분산이 매우 큰 경우 절사해서 비교
means = np.array(means)
medians = np.array(medians)
print(f"오염 비율 δ = {delta}")
print(f"표본 평균 — Var = {np.var(means):.4f}, "
f"90% 범위: [{np.percentile(means,5):.2f}, {np.percentile(means,95):.2f}]")
print(f"표본 중앙값 — Var = {np.var(medians):.6f}, "
f"90% 범위: [{np.percentile(medians,5):.3f}, {np.percentile(medians,95):.3f}]")
simulate_contamination(delta=0.0) # 오염 없음
print()
simulate_contamination(delta=0.05) # 5% Cauchy 오염
def breakdown_demo():
"""붕괴점 시연: 하나의 이상치가 평균에 미치는 영향."""
base_data = [1.2, 2.3, 3.1, 2.8, 1.9, 2.5, 3.0, 2.1]
print(f"\n원본 데이터 평균: {np.mean(base_data):.3f}, 중앙값: {np.median(base_data):.3f}")
for outlier in [10, 100, 1000, 1e6]:
data_with_outlier = base_data + [outlier]
print(f"이상치 = {outlier:.0e} | "
f"평균 = {np.mean(data_with_outlier):.2f}, "
f"중앙값 = {np.median(data_with_outlier):.2f}")
breakdown_demo()11.2 Step 2: Huber M-추정량과 ARE 계산 (실무 활용)
import numpy as np
from scipy.stats import norm, logistic, laplace
from scipy.optimize import brentq
# Huber 추정량 직접 구현
def huber_estimator(data, k=1.5, max_iter=100, tol=1e-8):
"""
Huber M-추정량을 반복 가중 최소제곱(IRLS)으로 계산.
Score 방정식: Σ ψ(xᵢ - θ) = 0 을 반복법으로 풀기.
"""
theta = np.median(data) # 초기값: 중앙값
for _ in range(max_iter):
residuals = data - theta
# Huber 가중치: |r| <= k이면 1, 아니면 k/|r|
weights = np.where(np.abs(residuals) <= k, 1.0, k / np.abs(residuals))
theta_new = np.sum(weights * data) / np.sum(weights)
if abs(theta_new - theta) < tol:
break
theta = theta_new
return theta
# 예시 데이터 (8개 정규 + 3개 이상치)
x = np.array([-1.28, -0.96, -0.46, -0.44, -0.26, -0.21, -0.063, 0.39, 3, 6, 9])
print("Huber 추정량 (다양한 k 값):")
for k_val in [0.001, 1, 3, 5, 10, 1e10]:
h_est = huber_estimator(x, k=k_val)
print(f" k = {k_val:8.3f} → Huber = {h_est:.4f}")
print(f" 중앙값 (k→0) = {np.median(x):.4f}")
print(f" 평균 (k→∞) = {np.mean(x):.4f}")
# ARE 계산: ARE(중앙값, 평균) = 4σ² f(μ)²
def are_median_mean(distribution='normal'):
"""
세 가지 분포에서 ARE(중앙값, 평균) 계산.
분자: σ² (평균의 점근 분산)
분모: 1/(4f(μ)²) (중앙값의 점근 분산)
"""
if distribution == 'normal':
f_mu = norm.pdf(0) # 1/√(2π)
sigma2 = 1.0 # N(0,1)의 분산
elif distribution == 'logistic':
f_mu = logistic.pdf(0) # 1/4
sigma2 = np.pi**2 / 3 # 로지스틱 분산
elif distribution == 'double_exponential':
f_mu = laplace.pdf(0) # 1/2
sigma2 = 2.0 # 이중지수 분산
var_mean = sigma2 # 평균의 점근 분산
var_median = 1 / (4 * f_mu**2) # 중앙값의 점근 분산
are = var_mean / var_median
print(f"[{distribution:20s}] f(μ) = {f_mu:.4f}, σ² = {sigma2:.4f}, "
f"Var(median) = {var_median:.4f}, ARE = {are:.4f}")
print("\nARE(중앙값, 평균) 비교:")
are_median_mean('normal')
are_median_mean('logistic')
are_median_mean('double_exponential')
# 붕괴점 vs ARE trade-off 시각화 (출력 기반)
print("\n추정량 비교 요약:")
print(f"{'추정량':<20} {'붕괴점':>8} {'정규 ARE':>10} {'강건성':>10}")
print("-" * 52)
print(f"{'표본 평균':<20} {'0%':>8} {'1.000':>10} {'낮음':>10}")
print(f"{'Huber (k=1.5)':<20} {'~7.5%':>8} {'0.960':>10} {'중간':>10}")
print(f"{'표본 중앙값':<20} {'50%':>8} {'0.637':>10} {'높음':>10}")12 핵심 요약
강건성 계층 (붕괴점 기준):
표본 평균 표본 중앙값
(BD = 0%) ←────────────────→ (BD = 50%)
Huber M-추정량
(BD = 0 ~ 50%)
↑ k 조절
k → 0: 중앙값에 수렴 (강건, 덜 효율적)
k → ∞: 평균에 수렴 (효율적, 덜 강건)
k = 1.5: 정규에서 4% 효율 손실, 이상치 내성 보유| 핵심 결과 | 내용 |
|---|---|
| 표본 평균의 취약점 | 코시 오염 → 분산 \(\infty\), 붕괴점 = 0 |
| 중앙값 점근 분산 | \(1/(4f(\mu)^2)\) |
| ARE(중앙값, 평균), 정규 | \(2/\pi \approx 0.637\) |
| M-추정량 ARE 상한 | \(\leq 1\) (MLE 대비, 코시-슈바르츠) |
| Huber \(k=1.5\), 정규 | ARE \(= 0.96\) — 효율성을 거의 포기하지 않고 강건성 획득 |
13 관련 주제
선행 지식
- 점근적 계산과 비교 (Asymptotic Calculations and Comparisons) — Delta Method, ARE 개념
- 부트스트랩 표준오차 (Bootstrap Standard Errors) — 분포 무관 분산 추정
후속 주제
- M-추정량 (M-Estimators) — Huber/Bisquare 비교, 영향함수, IRLS
- 점근적 가설검정 (Asymptotic Hypothesis Testing) — Wald, Score, LRT
관련 개념
- 순서통계량 (Order Statistics) — 중앙값의 분포, 분위수 이론
- 점근적 점추정: 효율성 — 크래머-라오 하한, MLE 효율성
14 참고 문헌
- Casella, G. & Berger, R. L. (2002). Statistical Inference (2nd ed.). Duxbury. Ch.10.
- Huber, P. J. (1964). Robust estimation of a location parameter. Annals of Mathematical Statistics, 35(1), 73–101.
- Huber, P. J. (1981). Robust Statistics. Wiley.
- Hampel, F. R. (1974). The influence curve and its role in robust estimation. JASA, 69, 383–393.