로버스트 분산추정량¶

import numpy as np
from sklearn.covariance import LedoitWolf
rng = np.random.default_rng(0)
p, n, R = 30, 60, 500
Sigma = 0.5**np.abs(np.subtract.outer(np.arange(p), np.arange(p)))
err_s = err_lw = 0.0; shr = []
for _ in range(R):
    X = rng.multivariate_normal(np.zeros(p), Sigma, n)
    S = np.cov(X, rowvar=False)
    lw = LedoitWolf().fit(X)
    err_s += np.linalg.norm(S - Sigma, "fro")**2
    err_lw += np.linalg.norm(lw.covariance_ - Sigma, "fro")**2
    shr.append(lw.shrinkage_)
print(f"MSE: sample={err_s/R:.1f}, LW={err_lw/R:.1f}, shrinkage={np.mean(shr):.2f}")

출력:

MSE: sample=16.2, LW=8.9, shrinkage=0.45

예상 결과: LW의 평균제곱오차가 대략 40–60% 낮고, \(p/n = 0.5\)에서 축소 강도는 \(\sim 0.3\)이다. \(p/n\)이 클수록 축소가 강해진다. 고차원 영역(\(p \approx n\))에서는 표본공분산이 특이행렬이 되는데, LW는 역행렬 존재성과 더 나은 평균제곱오차를 함께 준다.

중앙값 = 7. 절대편차: \(|2-7|, |4-7|, |6-7|, |8-7|, |10-7|, |100-7| = 5, 3, 1, 1, 3, 93\).

정렬하면 1, 1, 3, 3, 5, 93. 중앙값(MAD) = (3 + 3)/2 = 3.

축척된 MAD = \(1.4826 \cdot 3 \approx 4.45\).

표본표준편차 = \(\sqrt{(\sum(x_i - \bar x)^2)/(n-1)}\), \(\bar x = 21.67\). \(\sum(x_i - \bar x)^2 \approx 7403\). \(\mathrm{SD} \approx 38.5\).

표준편차는 이상점에 지배되는 반면, MAD는 자료 대부분의 전형적인 척도를 포착한다.

표본분산: 붕괴점 0%(오염된 관측값 하나가, 특히 제곱을 통해, 값을 얼마든지 옮긴다).

MAD: 붕괴점 50%. 중앙값으로부터의 절대편차의 중앙값 — 두 겹의 중앙값이 모두 로버스트하다.

IQR: 붕괴점 25%. 한쪽 꼬리의 25%를 오염시키면 사분위수 하나를 옮길 수 있다.

이상점이 잦은 자료에서는 MAD가 가장 로버스트하고 표본분산이 가장 취약하다.

\(n = 7\)에서 20% 절사: \(k = 1\). 가장 작은 값과 가장 큰 값을 절사하면 \(\{3, 5, 7, 9, 11\}\)이 남는다.

절사평균: 7 (앞의 연습문제).

절사분산: \(\frac{1}{n-2k}\sum_{i=k+1}^{n-k}(X_{(i)} - \bar X_{\text{trim}})^2 = (1/5)[(3-7)^2 + (5-7)^2 + (7-7)^2 + (9-7)^2 + (11-7)^2] = 40/5 = 8\).

원래 자료의 표본분산은 100에 지배된다. 절사분산은 전형적인 관측값들의 산포를 포착한다.

축척: 정규 자료에서 절사분산은 \(\sigma^2\)의 편향추정량이며, (MAD의 \(1.4826\)에 해당하는) 축척 상수가 존재한다.

표본공분산: \(S = (1/(n-1))\sum (X_i - \bar X)(X_i - \bar X)^\top\). 제곱편차가 어느 차원의 이상점이든 증폭한다.

대안:

• 최소공분산행렬식(MCD): 표본공분산의 행렬식이 가장 작은 \(\lceil n/2 \rceil\)개 관측값의 부분집합을 찾는다. 로버스트하지만 계산 비용이 크다.
• 최소부피타원체(MVE): \(\lceil n/2 \rceil\)개 관측값을 담는 가장 작은 타원체를 찾는다. 매우 로버스트하다.
• Tukey의 bisquare M-추정량: 중심으로부터의 Mahalanobis 거리에 따라 관측값의 가중치를 낮춘다.
• Ledoit-Wolf: 구조화된 목표(대각행렬) 쪽으로 축소한다. 이상점에는 로버스트하지 않지만 고차원 잡음을 다룬다.

다변량 이상점 탐지(금융 포트폴리오, 다변량 품질관리)에서는 MCD가 표준이다. sklearn.covariance.MinCovDet에 구현되어 있다.

로버스트 공분산을 쓸 때:

• 이상점 탐지: (로버스트한 \(\hat{\boldsymbol\Sigma}\)로 계산한 Mahalanobis 거리를 통해) 자료 대부분에서 멀리 떨어진 점을 찾고 싶을 때.
• 분류를 위한 Mahalanobis 거리: QDA, LDA, k-NN에 쓰는 공분산추정값을 이상점이 좌우해서는 안 될 때.
• 오염된 자료의 PCA: 로버스트 공분산이 로버스트한 주성분을 준다.
• 다변량 품질관리: 이상 패턴 탐지(Hotelling의 \(T^2\)).

표본공분산을 쓸 때:

• 자료가 깨끗할 때(통제된 실험, 모의실험).
• 축소추정량이 선호되는 고차원 / 저잡음 영역.
• 계산 효율이 필요할 때.

고차원에서의 주의: \(p > n\)이면 이상점과 무관하게 표본공분산이 특이행렬이 된다. 축소(Ledoit-Wolf)나 희소 방법(그래프 라소)이 필요하며, 로버스트 성분과 축소 성분이 둘 다 필요할 수도 있다.

MAD의 약점.

1. 효율이 낮다. 정규분포에서 점근효율이 37%에 그친다. 붕괴점 50%의 대가치고도 큰 손실이다.
2. 중심을 먼저 정해야 한다. \(\operatorname{median}|x_i-\tilde x|\)로 중앙값에 기대므로, 분포가 비대칭이면 편향이 생긴다.
3. 대칭을 암묵적으로 가정한다. 양쪽 꼬리가 다르면 한쪽만 반영된다.

\(S_n\) 추정량.

모든 쌍의 거리를 쓰므로 중심을 정할 필요가 없다.

• 붕괴점 50%, 정규 효율 58%.
• 비대칭분포에도 자연스럽게 적용된다.

\(Q_n\) 추정량.

쌍 거리의 약 1/4 분위수를 쓴다.

• 붕괴점 50%, 정규 효율 82%.
• 세 방법 중 효율이 가장 높다.

비교.

권고. \(Q_n\)이 현재 가장 좋은 절충으로 평가된다. 붕괴점 50%를 유지하면서 효율 82%를 낸다. 다만 MAD보다 계산이 복잡하고 구현이 덜 보편적이다.

MAD는 계산이 간단하고 널리 알려져 있다는 점에서 여전히 기본 진단 도구로 쓸 만하다. 이상치를 찾아내는 용도라면 효율이 중요하지 않기 때문이다. 척도 자체를 추정해 보고해야 한다면 \(Q_n\)이나 MM-추정 척도를 고려한다.

좌표별 접근의 한계. 각 변수마다 MAD를 계산해 대각행렬을 만들면 상관구조를 전혀 담지 못한다. 그런데 다변량 이상치의 상당수는 좌표별로는 평범한데 조합이 이상한 점이다.

예를 들어 키와 몸무게가 강하게 양의 상관인 자료에서 "키 190cm, 몸무게 50kg"인 사람은 각 변수만 보면 범위 안이지만 조합이 극도로 이례적이다. 좌표별 진단으로는 절대 잡히지 않는다.

MCD의 착상. \(n\)개 관측 중 \(h\)개(\(h\approx0.75n\))를 골라 그 부분집합의 공분산행렬식이 가장 작아지도록 한다.

그 부분집합의 평균과 공분산이 추정값이다.

왜 행렬식인가. 행렬식이 타원의 부피에 비례하므로, 가장 조밀하게 뭉친 \(h\)개를 찾는 것과 같다. 이상치는 어느 조밀한 부분집합에도 들어가지 못하므로 자동으로 배제된다.

성질.

• 붕괴점이 \((n-h+1)/n\)으로, \(h=\lceil(n+p+1)/2\rceil\)이면 약 50%다.
• 아핀 등변이다. 자료를 선형변환하면 추정값도 대응되게 변환된다. 좌표축의 선택에 의존하지 않는다는 뜻이며, 좌표별 접근에는 없는 성질이다.
• 일치성을 위해 보정 인수가 필요하다(부분집합만 쓰므로 분산이 과소추정된다).

계산. 모든 부분집합을 보는 것은 불가능하므로 FAST-MCD 알고리즘을 쓴다. 무작위 초기 부분집합에서 시작해 "마할라노비스 거리가 작은 \(h\)개로 갱신"을 반복하면 행렬식이 단조감소한다. 여러 초기값에서 돌려 최선을 고른다.

쓰임. 앞서 본 로버스트 마할라노비스 거리의 기반이며, 로버스트 주성분분석, 로버스트 판별분석, 회귀의 지렛값 진단에 쓰인다. sklearn.covariance.MinCovDet이 구현되어 있다.

문제. \(p/n = 0.5\)처럼 차원이 관측 수에 가까우면 표본공분산행렬 \(S\)가 심하게 왜곡된다.

구체적으로.

1. 고윳값이 퍼진다. 마르첸코-파스투르 법칙에 따라, 참 공분산이 단위행렬이어도 표본 고윳값이

$$ \left(1-\sqrt{p/n}\right)^2 \quad\text{부터}\quad \left(1+\sqrt{p/n}\right)^2 $$

까지 퍼진다. \(p/n=0.5\)면 \(0.086\)에서 \(2.91\)까지다. 참값이 모두 1인데 최대와 최소가 34배 차이 난다.

1. 가장 큰 고윳값은 과대, 가장 작은 것은 과소추정된다. 그런데 포트폴리오 최적화는 \(S^{-1}\)을 쓰므로 작은 고윳값이 증폭된다.

2. \(p>n\)이면 \(S\)가 특이해 역행렬이 아예 없다.

결과. 최소분산 포트폴리오 \(w\propto S^{-1}\mathbf{1}\)이 극단적인 가중치를 낸다. 어떤 자산에 \(+300\%\), 다른 자산에 \(-250\%\) 같은 값이 나오고, 표본 밖 성능이 형편없다. "오차 극대화 절차"라는 별명이 붙은 이유다.

르두아-울프의 해법. 구조화된 목표 \(F\)(예: 모든 상관이 같은 행렬, 또는 \(\bar\lambda I\))를 향해 축소한다.

그리고 최적 \(\delta\)를 자료에서 해석적으로 추정한다. 기대 제곱오차를 최소로 하는 값이며, \(p/n\)이 클수록 \(\delta\)가 커진다.

효과. 고윳값의 퍼짐이 줄어들어 역행렬이 안정되고, 포트폴리오 가중치가 합리적인 범위에 들어온다. 표본 밖 샤프비율이 크게 개선된다는 것이 실증적으로 확인되어 있다.

대안들. 요인모형(주성분이나 파마-프렌치 요인으로 공분산을 구조화), 비선형 축소(고윳값마다 다르게 축소), 가중치에 직접 제약(공매도 금지)도 비슷한 효과를 낸다. 공매도 금지 제약이 사실상 축소와 같은 역할을 한다는 결과가 알려져 있다.

왜 필요한가. 로버스트 척도는 그 자체로 \(\sigma\)를 추정하지 않는다. MAD는 "중앙값으로부터의 중앙 절대편차"이고, 이것이 \(\sigma\)와 같지 않다. 기준 분포(보통 정규)에서 \(\sigma\)와 같아지도록 상수를 곱해 준다.

계산 방법. 기준 분포에서 그 척도의 참값을 구하고 역수를 취한다.

MAD. \(X\sim N(\mu,\sigma^2)\)에서

이므로

을 곱한다.

IQR.

이므로 \(\text{IQR}/1.349\)가 \(\sigma\)의 추정값이다.

\(Q_n\)과 \(S_n\). 같은 방식으로 정규분포에서의 상수를 구한다. \(Q_n\)은 2.2219, \(S_n\)은 1.1926이다.

중요한 단서 세 가지.

1. 보정은 기준 분포에 묶여 있다. 1.4826은 정규분포에서만 \(\sigma\)의 일치추정량을 준다. 자료가 \(t_5\)라면 이 상수가 \(\sigma\)를 과소추정한다. 그러나 대개 이것이 문제가 아닌데, 로버스트 척도를 쓰는 목적이 "오염되지 않았다면 \(\sigma\)가 얼마였을까"를 추정하는 것이기 때문이다.

2. 소표본 보정이 따로 필요하다. 위 상수들은 점근적인 값이다. \(n\)이 작으면 추가 보정 인수가 필요하며, 표로 정리되어 있다. \(n=10\)에서 MAD의 소표본 보정이 10%를 넘는다.

3. 보정하지 않은 값을 보고하면 안 된다. MAD 자체는 \(\sigma\)의 0.67배이므로, 보정 없이 "표준편차"라고 부르면 30% 과소보고다.

실무 권고. 소프트웨어가 보정을 넣었는지 확인한다. scipy.stats.median_abs_deviation은 기본적으로 보정하지 않으며 scale='normal'을 주어야 1.4826이 곱해진다. statsmodels.robust.mad는 기본적으로 보정한다. 규약이 다르므로 반드시 문서를 확인해야 한다.