13회 · 최소제곱 · 노트

1 결합의 문제 (1749–1799)

역사 18세기 천문학과 측지학은 같은 벽에 부딪혔다. 이론은 몇 개의 상수로 정해지는 식을 주고, 관측은 그 식에 넣을 숫자를 준다. 관측이 미지수보다 많으면 방정식들이 서로 어긋난다. Euler(1749)는 목성–토성 문제의 75개 방정식 앞에서 멈췄다. Mayer(1750)는 방정식을 묶어 더했고, Boscovich(1757)는 보정의 절댓값 합을 최소로 했으며, Laplace(1787)는 Mayer와 비슷한 결합으로 토성의 “큰 부등”을 풀었다. 이 반세기 동안 “관측을 결합한다”는 생각 자체가 자리 잡았다. Stigler는 이것을 통계학의 첫째 기둥 Aggregation의 탄생으로 본다.

1.1 선형모형

정의 1.1 (선형모형) $n$개의 관측이

$$\cdat{\y}=X\cpar{\bbeta}+\beps,\qquad E\beps=\bzero,\qquad \Cov\beps=\cpar{\sigma^2}I_n$$

을 만족한다고 하자. $X$는 $n\times p$의 상수 행렬이고 $\rank X=p\le n$이다(완전 계수). $\bbeta\in\R^p$와 $\sigma^2>0$은 모르는 모수다. 오차에 대한 두 조건(평균 0, 서로 무상관이고 분산이 같음)을 Gauss–Markov 가정이라 한다. 정규분포는 가정하지 않는다.

완전 계수는 “미지수가 실제로 $p$개”라는 뜻이다. 어떤 열이 다른 열들의 일차결합이면 그 열에 해당하는 모수는 자료로 가려낼 수 없다. $p=1$, $X=\one$이면 $y_i=\mu+\varepsilon_i$, 곧 한 양을 $n$번 잰 모형이다. 1–12회의 “iid 표본”은 이 모형의 가장 단순한 경우다.

예제 1.2 (Mayer의 27개 방정식) Mayer가 $i$번째 날 잰 Manilius의 위치에서 계산한 값을 $y_i$(도)라 하면, 칭동 이론은

$$y_i=-\beta+\alpha\,x_{i2}+(\alpha\sin\theta)\,x_{i3}+\varepsilon_i\qquad(i=1,\dots,27)$$

를 준다. $\beta$는 달 적도에서 잰 Manilius의 위도이고, $\alpha$와 $\theta$는 달 적도면의 기울기를 정하는 두 각이다. $x_{i2},x_{i3}$은 그날 달의 위치에서 계산되는 삼각함수 값이다. 처음 세 날의 값은 다음과 같다.

$$\begin{array}{cccc} i & y_i & x_{i2} & x_{i3}\\\hline 1 & -13.167 & 0.8836 & -0.4682\\ 2 & -13.133 & 0.9996 & -0.0282\\ 3 & -13.200 & 0.9899 & \phantom{-}0.1421 \end{array}$$

$\theta$가 사인 속에 있어서 모형은 $(\beta,\alpha,\theta)$에 대해 선형이 아니다. 그러나 $\gamma=\alpha\sin\theta$로 바꿔 쓰면 $(\beta,\alpha,\gamma)$의 선형모형이 된다. $X$의 첫 열은 모두 $-1$, 둘째·셋째 열은 $x_{i2},x_{i3}$이다. $\gamma$를 추정한 뒤 $\hat\theta=\arcsin(\hat\gamma/\hat\alpha)$로 돌아간다. 9회 정리 1.3(MLE의 불변성)과 같은 재모수화다.

정의 1.3 (선형불편추정량) $\btil=A\y$ ($A$는 $p\times n$ 상수 행렬) 꼴의 추정량을 선형추정량이라 한다. 모든 $\bbeta$에서 $E\btil=\bbeta$이면 선형불편추정량(LUE)이다.

명제 1.4 (선형불편의 조건) 선형추정량 $A\y$가 불편 $\iff$ $AX=I_p$.

증명 $E A\y=AX\bbeta$이다(5회 보조정리 3.1). 이것이 모든 $\bbeta$에서 $\bbeta$와 같으려면 $\bbeta$에 표준기저 $\mathbf e_1,\dots,\mathbf e_p$를 차례로 넣어 $AX$의 열이 $I_p$의 열과 같아야 한다. 거꾸로 $AX=I_p$이면 $AX\bbeta=\bbeta$이다.

1.2 Mayer의 묶음 평균

Mayer는 27개 방정식을 $x_{i2}$(곧 $\alpha$의 계수)의 크기에 따라 9개씩 세 묶음으로 나눴다. 묶음 I은 $x_{i2}$가 1에 가까운 날(1, 2, 3, 6, 9, 10, 11, 12, 27번), 묶음 II는 $x_{i2}$가 음수인 날(8, 18, 19, 21–26번), 묶음 III은 나머지다. 묶음마다 방정식을 더하면

$$\begin{align*} -9\beta+8.4987\,\alpha-0.7932\,\gamma&=-118.133,\\ -9\beta-6.1404\,\alpha+1.7443\,\gamma&=-140.283,\\ -9\beta+2.7977\,\alpha+7.9649\,\gamma&=-127.533 \end{align*}$$

의 세 방정식이 된다. 오차항은 더한 뒤 버린다.

예제 1.5 (묶음 평균은 선형불편추정량이다) 묶음을 나타내는 $27\times3$ 행렬 $G$($i$가 묶음 $j$에 속하면 $G_{ij}=1$, 아니면 0)를 쓰면 묶음의 합은 $G\T\y=G\T X\bbeta+G\T\beps$이고, Mayer의 해는

$$\btil_{\mathrm M}=(G\T X)^{-1}G\T\y$$

이다. $A=(G\T X)^{-1}G\T$로 두면 $AX=I_3$이므로 명제 1.4로 선형불편이다. 위의 세 방정식을 풀면 $\hat\beta=14.547^\circ=14^\circ32.8'$, $\hat\alpha=1.4958^\circ=1^\circ29.7'$, $\hat\gamma=-0.0996$이고 $\hat\theta=\arcsin(\hat\gamma/\hat\alpha)=-3^\circ49.1'$이다(연습 1.1). 아무 세 방정식을 골라 푸는 것도 $G$가 선택 행렬인 경우이므로 역시 선형불편이다. 27개에서 3개를 고르는 방법은 $\binom{27}{3}=2925$가지다.

모든 선형불편추정량은 평균적으로 옳다. 그렇다면 무엇으로 고르나? 답은 분산이다. Mayer는 묶어 더하면 정확도가 결합한 방정식의 수에 비례해 좋아진다고 썼다. 실제로는 표준편차가 그 수의 제곱근에 비례해 줄어든다(연습 1.2). 비례 상수가 틀렸어도, “더하면 좋아진다”는 직관은 옳았다. 3절의 Gauss–Markov 정리는 모든 선형불편추정량 가운데 분산이 가장 작은 것이 최소제곱임을 보인다. Mayer의 해는 그보다 분산이 크다는 것도 거기서 계산한다.

1.3 Boscovich의 절댓값 (1757)

Newton은 지구가 극 쪽으로 납작하다고 예측했다. 그렇다면 위도 1도에 해당하는 자오선 호의 길이는 극으로 갈수록 길어야 한다. 위도 $L$에서의 길이는 근사적으로

$$a(L)=a_0+a_1\sin^2L$$

이고, 납작한 정도(편평률)는 $f\approx a_1/(3a_0)$이다. 1750년대까지 다섯 곳에서 이 길이를 쟀다. 미지수 2개에 방정식 5개다.

장소	위도	$x=\sin^2L$	$y$ (toise)	Boscovich 해의 잔차
Quito	$0^\circ0'$	0	56751	0.0
희망봉	$33^\circ18'$	0.2987	57037	79.2
로마	$42^\circ59'$	0.4648	56979	$-93.8$
파리	$49^\circ23'$	0.5762	57074	$-75.9$
Lapland	$66^\circ19'$	0.8386	57422	90.5

Boscovich는 직선 $y=b_0+b_1x$에 두 조건을 걸었다. (i) 보정(잔차)의 합이 0이다. (ii) 그 조건 아래 보정의 절댓값 합이 가장 작다.

명제 1.6 (Boscovich의 해는 가중 중앙값이다) $x_i$들이 모두 같지는 않다고 하자. 조건 (i)을 만족하는 직선은 무게중심 $(\bar x,\bar y)$를 지난다. 그 위에서 절댓값 합은

$$\sum_{i=1}^n|y_i-b_0-b_1x_i|=\sum_{i:\,x_i\ne\bar x}|x_i-\bar x|\,\bigl|s_i-b_1\bigr|+\sum_{i:\,x_i=\bar x}|y_i-\bar y|,\qquad s_i=\frac{y_i-\bar y}{x_i-\bar x}$$

이고, 이를 최소로 하는 $b_1$은 기울기 $s_i$들의 가중 중앙값(무게 $w_i=|x_i-\bar x|$)이다. 곧 $s_i$를 크기 순으로 늘어놓고 무게를 누적했을 때 누적 무게가 처음으로 전체의 절반 이상이 되는 $s_i$다.

증명 (i) $\sum(y_i-b_0-b_1x_i)=0$은 $b_0=\bar y-b_1\bar x$와 같다. 대입하면 $y_i-b_0-b_1x_i=(y_i-\bar y)-b_1(x_i-\bar x)$이고, $x_i\ne\bar x$인 항은 $|x_i-\bar x|\,|s_i-b_1|$로 쓸 수 있다. 나머지 항은 $b_1$과 무관하다.

$g(b)=\sum w_i|s_i-b|$는 볼록한 조각별 일차함수이고, $s_i$들이 아닌 점에서의 기울기는 $g'(b)=\sum_{s_i<b}w_i-\sum_{s_i>b}w_i$이다. 이 값은 $b$가 커지면서 $-\sum w_i$에서 $+\sum w_i$까지 늘어나고, 처음으로 0 이상이 되는 곳이 최소점이다. 그 점이 누적 무게가 절반을 넘는 $s_i$, 곧 가중 중앙값이다.

예제 1.7 (다섯 개의 호) 무게중심은 $(\bar x,\allowbreak \bar y)=(0.4357,\allowbreak \ 57052.6)$이다. 다섯 기울기를 크기 순으로 쓰면 로마 $-2525.7$ (무게 0.029), 희망봉 113.9 (0.137), 파리 152.3 (0.141), Quito 692.3 (0.436), Lapland 916.8 (0.403)이다. 누적 무게의 비율은 0.025, 0.145, 0.268, 0.648, 1이므로 가중 중앙값은 Quito의 기울기 692.3이다. Boscovich의 직선은 $y=56751+692.3\,x$이고 편평률은 $1/246$이다. 최소제곱 직선은 $y=56737.4+723.4\,x$ (편평률 $1/235$)이고, 조건 (i) 없이 절댓값 합만 최소로 하면 Quito와 Lapland를 지나는 $y=56751+800.1\,x$ ($1/213$)이다. 오늘의 값은 약 $1/298$이다. 다섯 측정 자체에 계통오차가 있었다(그림 1.1).

그림 1.1 Boscovich의 다섯 호. 왼쪽: 최소제곱(파랑), Boscovich(주황), 제약 없는 최소절댓값(보라)의 세 직선. 오른쪽: 무게중심을 지나는 직선들의 절댓값 합 $g(b_1)$. 꺾이는 점이 다섯 기울기 $s_i$이고, 최소점은 가중 중앙값인 Quito의 기울기 692.3이다.

가정 주의 표의 희망봉 값 0.2987은 Boscovich가 위도를 $33^\circ8'$로 잘못 넣어 계산한 것이다. $33^\circ18'$이면 0.3014다. 이 오류는 결론을 거의 바꾸지 않는다(최소제곱 편평률 $1/235.3\to1/234.9$). 이 노트는 역사적 계산을 재현하려고 원래 값을 쓴다.

조건 (ii)만 쓰는 최소절댓값(LAD) 추정은 1799년 Laplace가 알고리즘으로 정리했고, 1978년 Koenker와 Bassett가 분위수 회귀로 일반화했다. 다음 명제는 그 해가 왜 늘 몇 개의 관측을 정확히 지나는지 설명한다.

명제 1.8 (최소절댓값 해는 관측 $p$개를 지난다) $\rank X=p$이면 $\sum_i|y_i-\x_i\T\bb|$를 최소로 하는 $\bb$ 가운데, $\x_i$들이 일차독립인 관측 $p$개에서 잔차가 0인 것이 있다.

증명 $p$개의 초평면 $\{\bb:\x_i\T\bb=y_i\}$로 $\R^p$를 나누면, 각 조각 안에서 잔차의 부호가 고정되므로 목적함수는 일차함수다. $\x_i$들이 $\R^p$를 생성하므로 각 조각은 직선을 포함하지 않는 다면체이고, 꼭짓점은 일차독립인 초평면 $p$개의 교점이다. 목적함수는 0 이상이므로 각 조각에서 아래로 유계이고, 아래로 유계인 일차함수는 그런 다면체 위에서 꼭짓점에서 최솟값을 갖는다(선형계획법의 기본 정리). 전체 최솟값은 어느 조각의 최솟값이므로 어느 꼭짓점에서 달성된다.

다섯 호에서 제약 없는 해가 Quito와 Lapland를 지나는 것, Boscovich의 해가 무게중심과 Quito를 지나는 것이 이 명제의 모습이다. Mayer의 자료에서 최소절댓값 해는 $\hat\beta=14^\circ34.1'$, $\hat\alpha=1.515^\circ$이고 세 방정식을 정확히 지난다(6절).

연습문제

연습 1.1 (Mayer의 해 ★) 본문의 묶음 방정식 세 개를 풀어 예제 1.5의 $\hat\beta,\allowbreak \hat\alpha,\allowbreak \hat\gamma,\allowbreak \hat\theta$를 확인하라.

풀이

첫 식에서 둘째 식을 빼면 $14.6391\,\alpha-2.5375\,\gamma=22.150$, 셋째 식에서 둘째 식을 빼면 $8.9381\,\alpha+6.2206\,\gamma=12.750$이다. 이 둘을 풀면 $\hat\alpha=1.4958$, $\hat\gamma=-0.0996$이고, 첫 식에 넣으면 $-9\hat\beta=-118.133-8.4987\hat\alpha+0.7932\hat\gamma$에서 $\hat\beta=14.547$이다. $\hat\theta=\arcsin(-0.0996/1.4958)=-3.82^\circ=-3^\circ49.1'$.

연습 1.2 (정확도는 제곱근으로 ★) $\varepsilon_1,\allowbreak \dots,\allowbreak \varepsilon_m$이 서로 무상관이고 분산이 $\sigma^2$이면 평균 $\bar\varepsilon$의 표준편차가 $\sigma/\sqrt m$임을 보여라. Mayer의 묶음 하나($m=9$)를 더해 9로 나누면 표준편차는 몇 분의 1이 되는가?

풀이

무상관이므로 $\Var\bar\varepsilon=\tfrac1{m^2}\sum\Var\varepsilon_i=\sigma^2/m$이고 표준편차는 $\sigma/\sqrt m$이다. $m=9$이면 3분의 1이다. Mayer가 생각한 9분의 1이 아니다.

연습 1.3 (위치모형의 세 방법 ★) $y_i=\mu+\varepsilon_i$ ($X=\one$)에서 (a) 최소제곱 해는 $\bar y$, (b) 최소절댓값 해는 표본 중앙값, (c) Boscovich의 조건 (i)만으로 해가 $\bar y$로 정해짐을 보여라.

풀이

(a) 정리 2.2에서 $X=\one$이면 $n\hat\mu=\sum y_i$. (b) $\sum|y_i-\mu|$의 최소점은 중앙값이다(5회 연습 1.4). (c) 조건 (i)은 $\sum(y_i-b_0)=0$, 곧 $b_0=\bar y$이고, 미지수가 하나뿐이어서 조건 (ii)가 할 일이 남지 않는다.

연습 1.4 (선형계획법 ★★) 최소절댓값 문제가 선형계획법

$$\min\ \sum_{i=1}^n(u_i+v_i)\quad\text{s.t.}\quad X\bb+\bu-\bv=\y,\ \ \bu,\bv\ge\bzero$$

과 같음을 보여라. 최적해에서 각 $i$마다 $u_iv_i=0$임을 보여라.

풀이

$\bb$를 고정하면 $u_i-v_i=y_i-\x_i\T\bb$를 만족하는 $u_i,v_i\ge0$ 가운데 $u_i+v_i$의 최소는 $|y_i-\x_i\T\bb|$이고 $u_i=\max(r_i,0)$, $v_i=\max(-r_i,0)$일 때 이루어진다. 따라서 두 문제의 최솟값과 최소점의 $\bb$가 같다. 최적해에서 $u_i,v_i>0$이면 둘에서 $\min(u_i,v_i)$를 빼도 제약이 유지되고 목적함수가 작아지므로 모순이다.

연습 1.5 (세 방정식만 고르기 ★★) Mayer의 방정식에서 세 개를 골라 정확히 풀어 얻는 추정량을 $\btil_S=X_S^{-1}\y_S$라 하자($X_S$는 고른 세 행). $AX=I_3$인 $A$를 써서 이것이 선형불편임을 보이고, $\Cov\btil_S=\sigma^2(X_S\T X_S)^{-1}$임을 보여라. 2925가지 가운데 $\det X_S$가 0에 가까운 것은 왜 나쁜가?

풀이

$E_S$를 고른 세 행을 뽑는 $3\times27$ 행렬이라 하면 $\y_S=E_S\y$, $X_S=E_SX$이고 $A=X_S^{-1}E_S$에 대해 $AX=X_S^{-1}X_S=I_3$이다. $\Cov\btil_S=\sigma^2AA\T=\sigma^2X_S^{-1}E_SE_S\T X_S^{-\mathsf T}=\sigma^2(X_S\T X_S)^{-1}$. $\det X_S$가 0에 가까우면 $(X_S\T X_S)^{-1}$이 커져 분산이 폭발한다. 2925가지 해의 $\hat\beta$가 $-426^\circ$에서 $104^\circ$까지 흩어진 것이 그 결과다(노트북).

이 절을 PDF로 보기 (3쪽)노트북 1절

이 회차에서 읽은 절 0개 · 회차 차례