ReflexKR

드론 관제실의 수치해석 기초: NumPy로 근사해 풀기

강의 개요

공식이 없으면 한 걸음씩: 수치해석이란

0.1 + 0.2가 0.3이 아닌 날: 컴퓨터 속 실수

근사의 대가: 테일러 급수와 절단 오차

구간을 반으로: 이분법과 가위치법

접선을 타고 내려가기: 뉴턴-랩슨과 할선법

가장 멀리 나는 속도: 1차원 최적화

배송 허브 세 곳의 물량 맞추기: 가우스 소거와 피벗

지도가 커지면: LU 분해, 반복법, 조건수

센서 기록에 선 긋기: 최소제곱 회귀

비어 있는 시각 채우기: 보간과 스플라인

얼마나 날았나: 수치 적분

처음 문제로 돌아가기: 미분방정식의 수치해

드론 관제실의 수치해석 기초: NumPy로 근사해 풀기

9교시

강의 시간 6:40–7:30

센서 기록에 선 긋기: 최소제곱 회귀

운항 일지

#09 · 배터리 기록

남은 비행 시간은 몇 분일까

문제

섬 노선을 나는 6호기의 배터리 기록이 2분마다 들어온다. 센서 잡음 때문에 숫자가 들쭉날쭉해서, 예비 배터리 20%에 언제 닿을지 바로 읽을 수 없다.

목표

기록 전체에 가장 잘 맞는 직선과 곡선을 최소제곱으로 구하고, 그 예측을 얼마나 믿을 수 있는지(R², 표준 오차, 외삽의 위험) 숫자로 말한다.

하늘상자 새벽 배송 관제실

입사 둘째 주, 관제실 화면에 6호기의 배터리 기록이 2분마다 한 줄씩 올라옵니다. 민재 팀장이 화면을 가리킵니다. "26분째인데 41.9%예요. 예비 20%는 남기고 내려야 하거든요. 몇 분 더 날 수 있어요? 섬 마을에 한 군데 더 들를지 지금 정해야 해요."

서윤은 마지막 두 기록으로 기울기를 잽니다. 24분에 47.6%, 26분에 41.9%이니 1분에 2.85%씩 줄었습니다. 남은 21.9%를 2.85로 나누면 7.7분입니다. 그런데 바로 앞 두 기록(22분 52.6%, 24분 47.6%)으로 재면 1분에 2.5%, 남은 시간은 8.8분입니다. "어느 두 점을 고르느냐에 따라 1분 넘게 달라지네요."

오 박사가 고개를 끄덕입니다. "센서 기록에는 잡음이 끼어 있어. 두 점만 믿으면 잡음 두 개를 믿는 셈이지. 점을 전부 쓰는 선을 그어. 대신 '가장 잘 맞는다'가 무슨 뜻인지부터 정해야 해."

이번 교시는 잡음 낀 기록에 모델을 맞추는 가장 기본 도구, 최소제곱 회귀(least-squares regression) 를 다룹니다. 선형대수 강의 13교시에서는 최소제곱을 사영이라는 그림으로 봤습니다. 여기서는 계산하는 쪽에서 봅니다. 공식과 NumPy 함수, 맞춘 결과를 평가하는 숫자, 그리고 계산이 망가지는 곳까지 다룹니다.

9-1. 기록 위에 선 하나: 잔차와 제곱합 (7분)

먼저 기록을 만듭니다. 실제 기록 대신, 숨은 모델에 센서 잡음(표준편차 0.6%)을 더해 만든 숫자를 씁니다. 시드를 고정했으니 누구나 똑같은 기록을 얻습니다. 숨은 모델은 우리만 아는 정답이고, 관제실은 모른다고 칩니다.

>>> rng = np.random.default_rng(7)
>>> t = np.arange(0, 27, 2)          # 비행 시간 (분), 2분마다
>>> y = np.round(100 - 1.9*t - 0.012*t**2 + rng.normal(0, 0.6, t.size), 1)
>>> t
array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22, 24, 26])
>>> y
array([100. ,  96.3,  92. ,  87.6,  83.8,  79.2,  75.5,  71.9,  66.2,
        61.5,  57.5,  52.6,  47.6,  41.9])

직선 $\hat{y} = a + b\,t$ 를 하나 그었다고 합시다. $\hat{y}$ (와이 햇)은 "모델이 예측한 값"이라는 뜻입니다. 기록 $y_i$ 와 예측 $\hat{y}_i$ 의 차이를 잔차(residual) 라고 합니다.

r_i = y_i - \hat{y}_i = y_i - (a + b\,t_i)

잔차는 그래프에서 점과 직선 사이의 세로 거리입니다. 가로나 수직 거리가 아닌 이유는, 시간 $t$ 는 정확하고 잡음은 배터리 값 $y$ 에만 있다고 보기 때문입니다.

서윤이 눈대중으로 그은 직선 $\hat{y} = 100 - 2.2\,t$ 의 잔차를 봅시다. 잔차를 한 숫자로 모으는 방법은 세 가지가 떠오릅니다.

>>> r = y - (100 - 2.2 * t)
>>> np.round(r, 1)
array([ 0. ,  0.7,  0.8,  0.8,  1.4,  1.2,  1.9,  2.7,  1.4,  1.1,  1.5,
        1. ,  0.4, -0.9])
>>> round(r.sum(), 2), round(np.abs(r).sum(), 2), round((r**2).sum(), 2)
(np.float64(14.0), np.float64(15.8), np.float64(23.46))

후보	식	문제
잔차의 합	$\sum r_i$	+와 −가 서로 지운다. 평균 높이의 가로줄 $\hat{y} = \bar{y}$ 도 합이 0이다
절댓값의 합	$\sum \lvert r_i \rvert$	꺾인 점이 있어 미분으로 풀 수 없고, 답이 여러 개일 때도 있다
제곱의 합	$\sum r_i^2$	매끄러워서 미분 = 0으로 공식이 나온다. 큰 잔차를 더 무겁게 벌한다

그래서 잔차 제곱합(sum of squared residuals) 을 가장 작게 만드는 직선을 고릅니다. 이것이 최소제곱입니다. 잡음이 정규분포를 따른다면 최소제곱 직선이 "그 기록이 나왔을 가능성이 가장 큰 직선"이기도 합니다. 대신 제곱 때문에 크게 벗어난 점(이상치) 하나에 끌려가기 쉽습니다.

S(a, b) = \sum_{i=1}^{n} \bigl(y_i - a - b\,t_i\bigr)^2 \quad\text{을 가장 작게}

9-2. 직선 회귀: 공식, 정규방정식, NumPy (10분)

$S(a, b)$ 는 $a$ 와 $b$ 에 대한 아래로 볼록한 2차식입니다. 가장 낮은 곳에서는 두 방향의 기울기가 모두 0입니다. $a$ 로 미분한 것과 $b$ 로 미분한 것을 0으로 놓고 정리하면 연립방정식 두 개가 나옵니다.

\begin{aligned} n\,a + \Bigl(\sum t_i\Bigr) b &= \sum y_i \\ \Bigl(\sum t_i\Bigr) a + \Bigl(\sum t_i^2\Bigr) b &= \sum t_i y_i \end{aligned}

이것을 풀면 교과서의 공식이 됩니다. $\bar{t}, \bar{y}$ 는 평균입니다.

b = \frac{n \sum t_i y_i - \sum t_i \sum y_i}{n \sum t_i^2 - \bigl(\sum t_i\bigr)^2}, \qquad a = \bar{y} - b\,\bar{t}

손으로 한 번 넣어 봅시다. 기록에서 $n = 14$ , $\sum t_i = 182$ , $\sum y_i = 1013.6$ , $\sum t_i^2 = 3276$ , $\sum t_i y_i = 11166.4$ 입니다.

b = \frac{14 \times 11166.4 - 182 \times 1013.6}{14 \times 3276 - 182^2} = \frac{-28145.6}{12740} \approx -2.2092, \qquad a = 72.4 + 2.2092 \times 13 \approx 101.12

1분에 약 2.21%씩 줍니다. 400 Wh 배터리의 2.21%는 8.8 Wh이니, 평균 530 W쯤으로 날고 있다는 뜻입니다.

같은 연립방정식을 행렬로 적으면 쓰임이 넓어집니다. 열이 $[1, t_i]$ 인 행렬 $A$ 를 만들면 "모든 점을 지나는 직선" $A\mathbf{c} = \mathbf{y}$ 는 식이 14개, 모르는 수가 2개라 풀리지 않습니다. 양변에 $A^{\mathsf T}$ 를 곱한 2×2 연립방정식이 바로 위의 두 식입니다. 이것을 정규방정식(normal equations) 이라고 합니다.

A^{\mathsf T} A\,\mathbf{c} = A^{\mathsf T}\mathbf{y}, \qquad A = \begin{bmatrix} 1 & t_1 \\ 1 & t_2 \\ \vdots & \vdots \\ 1 & t_n \end{bmatrix},\quad \mathbf{c} = \begin{bmatrix} a \\ b \end{bmatrix}

>>> n = t.size
>>> A = np.column_stack([np.ones(n), t])
>>> A.T @ A
array([[  14.,  182.],
       [ 182., 3276.]])
>>> np.round(A.T @ y, 1)
array([ 1013.6, 11166.4])
>>> c = np.linalg.solve(A.T @ A, A.T @ y)
>>> np.round(c, 4)
array([101.12  ,  -2.2092])

$A^{\mathsf T}A$ 의 칸이 손 계산에 쓴 합들( $n$ , $\sum t_i$ , $\sum t_i^2$ )과 똑같습니다. 실무에서는 이 계산을 두 함수로 합니다.

>>> np.round(np.polyfit(t, y, 1), 4)
array([ -2.2092, 101.12  ])
>>> c2, ssr, rank, sv = np.linalg.lstsq(A, y)
>>> np.round(c2, 4), np.round(ssr, 3)
(array([101.12  ,  -2.2092]), array([9.382]))

np.polyfit(t, y, 차수)는 다항식을 맞추는 지름길이고, 계수를 높은 차수부터 돌려줍니다(기울기가 먼저). np.linalg.lstsq(A, y)는 어떤 $A$ 든 받는 일반 도구로, 잔차 제곱합과 계수(rank), 특이값도 함께 돌려줍니다. 둘 다 안에서는 정규방정식을 만들지 않고 $A$ 를 직접 분해합니다. 이유는 9-6에서 봅니다.

이제 민재 팀장의 질문에 답할 수 있습니다. 직선이 20%에 닿는 시각은 $a + b\,t = 20$ 을 풀면 됩니다.

>>> round((20 - c[0]) / c[1], 1)
np.float64(36.7)

36.7분, 지금(26분)부터 약 10.7분입니다. 두 점으로 잰 7.7분, 8.8분과 꽤 다릅니다. 그런데 이 직선을 믿어도 될까요?

9-3. 얼마나 잘 맞나: R²와 표준 오차 (6분)

맞춘 결과를 평가하는 숫자 두 개를 계산합니다. 먼저 기록이 평균 주위로 얼마나 퍼져 있는지 재는 총제곱합 $S_t$ 와, 직선이 설명하지 못하고 남긴 잔차 제곱합 $S_r$ 입니다.

S_t = \sum (y_i - \bar{y})^2, \qquad S_r = \sum (y_i - \hat{y}_i)^2, \qquad R^2 = 1 - \frac{S_r}{S_t}, \qquad s_e = \sqrt{\frac{S_r}{n - 2}}

결정계수(coefficient of determination) $R^2$ 는 "평균 하나로 때울 때보다 퍼짐을 몇 %나 줄였나"입니다. 1이면 모든 점을 지나고, 0이면 평균선보다 나을 것이 없습니다. 표준 오차(standard error) $s_e$ 는 잔차의 전형적인 크기로, 단위가 $y$ 와 같습니다(여기서는 %). $n$ 대신 $n-2$ 로 나누는 것은 계수 두 개를 기록에서 뽑아 썼기 때문입니다(자유도).

>>> yhat = A @ c
>>> SSR = np.sum((y - yhat)**2)
>>> SST = np.sum((y - y.mean())**2)
>>> round(SSR, 3), round(SST, 2)
(np.float64(9.382), np.float64(4450.82))
>>> round(1 - SSR / SST, 4)
np.float64(0.9979)
>>> round(np.sqrt(SSR / (n - 2)), 3)
np.float64(0.884)
>>> np.round(y - yhat, 1)
array([-1.1, -0.4, -0.3, -0.3,  0.4,  0.2,  0.9,  1.7,  0.4,  0.1,  0.6,
        0.1, -0.5, -1.8])

$R^2 = 0.998$ 이면 거의 완벽해 보입니다. 하지만 마지막 줄의 잔차 부호를 보세요. 양 끝은 −, 가운데는 +입니다. 잡음이라면 +와 −가 아무렇게나 섞여야 하는데, 여기서는 무늬가 있습니다. 실제 기록이 직선보다 살짝 아래로 휘어 있다는 신호입니다. 배터리가 줄수록 전압이 내려가 같은 전력에 전류가 커지므로, 방전이 조금씩 빨라지는 것이 자연스럽습니다.

표준 오차 0.884%도 단서입니다. 센서 잡음(0.6%)보다 큽니다. 직선이 잡음만이 아니라 휘어짐까지 잔차로 남겼다는 뜻입니다. $R^2$ 는 "얼마나 많이 설명했나"만 말하고 "모델 모양이 맞나"는 말하지 않습니다. 잔차를 늘 눈으로 보세요.

9-4. 곡선도 최소제곱으로: 다항 회귀와 선형화 (9분)

휘어진 기록에는 2차식 $\hat{y} = c_0 + c_1 t + c_2 t^2$ 을 맞춥니다. 곡선이지만 모르는 계수에 대해서는 일차식입니다. 그래서 $A$ 에 $t^2$ 열 하나만 더하면 똑같은 최소제곱 문제가 됩니다. 이렇게 열이 $1, t, t^2, \dots, t^m$ 인 행렬을 반데르몽드 행렬(Vandermonde matrix) 이라고 합니다.

A = \begin{bmatrix} 1 & t_1 & t_1^2 \\ \vdots & \vdots & \vdots \\ 1 & t_n & t_n^2 \end{bmatrix}, \qquad s_e = \sqrt{\frac{S_r}{n - (m + 1)}} \quad (m\text{차식은 계수가 } m+1\text{개})

>>> c_2 = np.polyfit(t, y, 2)
>>> c_2
array([-1.24484890e-02, -1.88557005e+00,  9.98253571e+01])
>>> for d in [1, 2, 3]:
...     cf = np.polyfit(t, y, d)
...     ssr = np.sum((y - np.polyval(cf, t))**2)
...     print(d, round(ssr, 3), round(1 - ssr / SST, 4), round(np.sqrt(ssr / (n - d - 1)), 3))
1 9.382 0.9979 0.884
2 2.162 0.9995 0.443
3 1.674 0.9996 0.409

np.polyval(계수, t)는 polyfit이 준 계수(높은 차수부터)로 값을 계산합니다. 2차식의 계수 $-0.01245$ 는 숨은 모델의 $-0.012$ 와 거의 같습니다. 1차 → 2차에서 표준 오차가 0.884 → 0.443으로 반이 되고, 2차 → 3차에서는 0.443 → 0.409로 별로 줄지 않습니다. 차수를 올려 표준 오차가 크게 줄면 그 휘어짐은 진짜이고, 조금만 줄면 잡음을 쫓기 시작한 것입니다. 2차가 알맞습니다.

2차식이 20%에 닿는 시각은 $c_2 t^2 + c_1 t + (c_0 - 20) = 0$ 의 근입니다. np.roots는 계수(높은 차수부터)를 받아 근을 모두 돌려줍니다.

>>> np.round(np.roots(c_2 - [0, 0, 20]), 1)
array([-186. ,   34.5])
>>> p2 = np.polynomial.Polynomial.fit(t, y, 2)
>>> round(p2(34.5), 2)
np.float64(19.96)

음수 근은 뜻이 없으니 34.5분, 지금부터 약 8.5분입니다. 직선의 10.7분보다 2분 넘게 짧습니다. 민재 팀장에게 2분은 마을 한 곳입니다. np.polynomial.Polynomial.fit은 같은 일을 하는 새 방식입니다. 결과를 함수처럼 부를 수 있고(p2(34.5)), 안에서 $t$ 를 $[-1, 1]$ 로 옮겨 맞추므로 차수가 높을 때 더 안전합니다(9-6).

선형화(linearization). 모델이 계수에 대해 일차식이 아니어도, 바꿔 적으면 일차식이 되는 경우가 있습니다. 착륙 뒤 배터리가 식는 기록이 그렇습니다. 바깥 기온과의 온도 차 $\Delta T$ 는 지수함수 $\Delta T = a\,e^{bt}$ 로 줄어듭니다( $b < 0$ ). 양변에 로그를 씌우면 $\ln \Delta T = \ln a + b\,t$ , 즉 $t$ 와 $\ln \Delta T$ 사이의 직선입니다.

>>> rng = np.random.default_rng(3)
>>> tm = np.arange(0, 31, 3)           # 착륙 뒤 시간 (분), 3분마다
>>> dT = np.round(18 * np.exp(-0.09 * tm) + rng.normal(0, 0.3, tm.size), 1)
>>> dT
array([18.6, 13. , 10.6,  7.8,  6. ,  4.6,  3. ,  2.6,  1.8,  2.6,  1.3])
>>> b_log, lna_log = np.polyfit(tm, np.log(dT), 1)
>>> round(np.exp(lna_log), 2), round(b_log, 4)
(np.float64(16.83), np.float64(-0.0847))

숨은 값은 $a = 18$ , $b = -0.09$ 인데 꽤 빗나갔습니다. 함정이 있습니다. 로그를 씌운 뒤의 최소제곱은 로그 공간의 잔차를 줄입니다. 로그는 작은 값의 차이를 크게 부풀립니다. 1.8과 2.6의 차이(0.8 °C)가 로그로는 $\ln 2.6 - \ln 1.8 \approx 0.37$ 로, 18.6과 13.0의 차이(5.6 °C)를 로그로 잰 0.36과 맞먹습니다. 잡음이 상대적으로 큰 꼬리 쪽 점들이 직선을 끌고 다닙니다.

로그 공간의 오차는 대략 (원래 오차) ÷ $\Delta T$ 입니다. 그래서 잔차마다 $\Delta T$ 를 곱해 원래 크기로 되돌리면 됩니다. np.polyfit의 w(가중치)가 바로 "잔차에 곱할 수"입니다.

>>> b_w, lna_w = np.polyfit(tm, np.log(dT), 1, w=dT)
>>> round(np.exp(lna_w), 2), round(b_w, 4)
(np.float64(18.06), np.float64(-0.0903))
>>> def ssr_exp(a, b):
...     return round(np.sum((dT - a * np.exp(b * tm))**2), 2)
>>> ssr_exp(np.exp(lna_log), b_log), ssr_exp(np.exp(lna_w), b_w)
(np.float64(4.84), np.float64(2.41))

가중치를 준 쪽이 숨은 값에 훨씬 가깝고, 원래 단위(°C)의 잔차 제곱합도 반으로 줄었습니다. 자주 쓰는 선형화를 모아 둡니다.

모델	바꿔 적기	`polyfit`에 넣을 것
지수 $y = a e^{bx}$	$\ln y = \ln a + b x$	`x`, `np.log(y)` (가중치 `w=y`)
거듭제곱 $y = a x^{b}$	$\ln y = \ln a + b \ln x$	`np.log(x)`, `np.log(y)`
포화 $y = \dfrac{a x}{b + x}$	$\dfrac{1}{y} = \dfrac{1}{a} + \dfrac{b}{a}\cdot\dfrac{1}{x}$	`1 / x`, `1 / y`

모두 $y > 0$ 이어야 하고, 바꿔 적은 공간에서 최소제곱을 한다는 점을 잊지 마세요.

9-5. 차수를 올리면 생기는 일: 과적합과 외삽 (8분)

차수를 올리면 잔차 제곱합은 절대 커지지 않습니다. 높은 차수의 다항식은 낮은 차수의 다항식을 모두 품고 있으니(높은 계수를 0으로 두면 된다), 최소제곱 답이 더 나빠질 수 없습니다. 그래서 잔차 제곱합이나 $R^2$ 만 보면 늘 "차수가 높을수록 좋다"는 답이 나옵니다. 함정입니다.

아래 실습대의 검은 점 14개가 9-1의 배터리 기록입니다. 남색 곡선이 최소제곱 다항식, 진홍 세로 선이 잔차입니다. 오른쪽 회색 띠는 기록이 없는 시간(외삽), 회색 점선은 예비 20%, 호박색 점은 곡선이 20%에 닿는 예측 시각입니다. 차수 슬라이더로 1차부터 9차까지 바꾸고, 점을 끌어(또는 점을 고르고 화살표 키로) 기록을 바꿔 보세요. '과적합 보기'를 켜면 2·8·14·20분의 점 4개가 청록으로 바뀌어 맞춤에서 빠지고, 곡선이 처음 보는 이 점들을 얼마나 잘 맞히는지(시험 점 RMSE) 따로 나옵니다. '차수 1 → 9 재생'은 차수를 올리며 곡선이 휘어 가는 모습을 이어서 보여 줍니다. 아래 NumPy 칸은 언제나 점 14개를 모두 써서 맞춘 결과입니다.

다항 회귀 곡선
잔차
기록 점 (끌 수 있음)
예비 20%

차수1잔차 보기과적합 보기 (점 4개를 시험용으로 빼기)

1차 다항식
잔차 제곱합: 9.382
결정계수 R²: 0.9979
표준 오차: 0.884 %

20%에 닿는 예측: 36.7분

차수가 낮아 곡선이 기록의 큰 흐름만 따라갑니다.

NumPy로 확인하기

실습대에서 값을 정할 때마다, 서버의 NumPy가 같은 숫자로 계산한 결과입니다.

>>> import numpy as np
>>> t = np.array([0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26])
>>> y = np.array([100, 96.3, 92, 87.6, 83.8, 79.2, 75.5, 71.9, 66.2, 61.5, 57.5, 52.6, 47.6, 41.9])
>>> c = np.polyfit(t, y, 1)
# 계수 (높은 차수부터)
>>> c
array([ -2.20923077, 101.12      ])
>>> res = y - np.polyval(c, t)
# 잔차 제곱합
>>> np.round(np.sum(res**2), 3)
np.float64(9.382)
# 결정계수 R²
>>> np.round(1 - np.sum(res**2) / np.sum((y - y.mean())**2), 4)
np.float64(0.9979)
# 기록 밖 34분의 예측 (외삽: 차수가 높을수록 믿기 어렵다)
>>> np.round(np.polyval(c, 34), 1)
np.float64(26.0)

차수를 1 → 2 → 3으로 올리며 R²와 20% 도달 예측이 어떻게 바뀌는지 보세요. 그다음 '과적합 보기'를 켜고 재생해 보세요. 9차가 되면 훈련 점은 모두 지나(RMSE 0) 시험 점에서는 크게 빗나가고, 회색 띠 안에서 곡선이 하늘로 솟습니다. 점 하나를 위로 끌어 이상치를 만들면 낮은 차수의 직선도 그 점 쪽으로 끌려갑니다.

실습대에서 본 것을 숫자로 확인합니다. 14개 중 4개(2·8·14·20분)를 시험용(held-out) 으로 떼어 두고, 나머지 10개로만 맞춘 뒤 두 묶음의 RMSE(잔차 제곱 평균의 제곱근)를 견줍니다. Polynomial.fit을 쓰는 것은 9차에서도 계산이 안정적이기 때문입니다.

>>> test = np.array([1, 4, 7, 10])
>>> train = np.setdiff1d(np.arange(n), test)
>>> for d in [1, 2, 3, 6, 9]:
...     p = np.polynomial.Polynomial.fit(t[train], y[train], d)
...     rmse_tr = np.sqrt(np.mean((y[train] - p(t[train]))**2))
...     rmse_te = np.sqrt(np.mean((y[test] - p(t[test]))**2))
...     print(d, round(rmse_tr, 3), round(rmse_te, 3))
1 0.728 1.096
2 0.285 0.663
3 0.237 0.618
6 0.126 0.614
9 0.0 1.409

훈련 RMSE는 차수와 함께 줄다가 9차에서 0이 됩니다. 점 10개를 계수 10개짜리 다항식이 모두 지나기 때문입니다. 시험 RMSE는 2차에서 크게 줄고, 그 뒤로는 제자리걸음을 하다가 9차에서 두 배로 뜁니다. 잡음까지 외운 곡선이 처음 보는 점에서 틀리는 것, 이것이 과적합(overfitting) 입니다. 시험 오차가 더 줄지 않는 가장 낮은 차수(여기서는 2차)를 고르는 것이 안전합니다.

과적합보다 더 무서운 것은 외삽(extrapolation), 기록 밖을 예측하는 일입니다. 민재 팀장의 질문(20%에 언제 닿나)은 처음부터 외삽입니다. 기록 14개를 모두 써서 맞춘 다항식으로 30분과 34분의 배터리를 예측해 봅니다.

>>> for d in [1, 2, 3, 6, 9]:
...     p = np.polynomial.Polynomial.fit(t, y, d)
...     print(d, round(p(30), 1), round(p(34), 1))
1 34.8 26.0
2 32.1 21.3
3 30.7 18.1
6 17.3 -58.8
9 131.8 1335.5

기록을 만든 숨은 모델로는 34분에 21.5%입니다. 2차가 가장 가깝고, 6차는 음수, 9차는 1335%를 내놓습니다. 기록 안에서는 9차가 가장 잘 맞았는데도 말입니다. 다항식은 기록 밖으로 나가면 가장 높은 차수 항이 모든 것을 정하기 때문에, 차수가 높을수록 기록 끝에서 빠르게 휩니다.

2차라고 안심할 수도 없습니다. 실제 리튬 배터리는 바닥 근처에서 전압이 갑자기 꺾여 떨어지는데, 26분까지의 기록에는 그 꺾임이 없습니다. 기록 안의 모양을 기록 밖으로 늘이는 것은 언제나 가정입니다. 그래서 관제실은 예측을 한 번 내고 끝내지 않고, 2분마다 새 기록이 들어올 때마다 다시 맞추고 여유를 둡니다.

오 박사의 한마디

R²가 0.999라고 기뻐하기 전에 두 가지를 물어. 잔차에 무늬가 있나? 그리고 지금 기록 밖을 예측하고 있나? 차수를 올려서 얻은 R²는 공짜가 아니야. 기록 밖에서 그 값을 치르지.

9-6. 계산이 망가지는 곳: 정규방정식과 조건수 (5분)

정규방정식 $A^{\mathsf T}A\,\mathbf{c} = A^{\mathsf T}\mathbf{y}$ 는 식은 깔끔하지만 계산에는 위험합니다. 8교시에서 배운 조건수(condition number)를 떠올려 봅시다. 조건수가 $10^k$ 이면 답에서 대략 $k$ 자리를 잃습니다. 그런데 $A^{\mathsf T}A$ 의 조건수는 $A$ 의 조건수의 제곱입니다. 잃는 자릿수도 두 배가 됩니다.

반데르몽드 행렬은 차수가 높아지면 조건수가 빠르게 커집니다. $t$ 를 $u = t/26$ 으로 $[0, 1]$ 에 옮겨도, $u^8$ 과 $u^9$ 처럼 높은 거듭제곱들은 모양이 서로 비슷해서(0 근처에서 납작하고 1 근처에서 솟는다) 열들이 거의 평행해지기 때문입니다.

>>> u = t / 26
>>> for d in [3, 6, 9, 11]:
...     V = np.vander(u, d + 1)
...     print(d, f"{np.linalg.cond(V):.0e}", f"{np.linalg.cond(V.T @ V):.0e}")
3 1e+02 1e+04
6 2e+04 3e+08
9 5e+06 3e+13
11 3e+08 2e+17

11차에서 $A^{\mathsf T}A$ 의 조건수가 $10^{17}$ , float64의 자릿수(약 16자리)를 다 잃는 크기입니다. 정규방정식과 np.linalg.lstsq(정규방정식을 만들지 않고 $A$ 를 직접 분해)로 같은 문제를 풀어 잔차 제곱합을 견줘 봅니다.

>>> for d in [9, 11, 13]:
...     V = np.vander(u, d + 1)
...     c_ne = np.linalg.solve(V.T @ V, V.T @ y)
...     c_ls = np.linalg.lstsq(V, y)[0]
...     print(d, round(np.sum((y - V @ c_ne)**2), 3), round(np.sum((y - V @ c_ls)**2), 3))
9 0.885 0.885
11 0.861 0.203
13 0.394 0.0

9차까지는 둘이 같습니다. 11차부터 정규방정식이 무너집니다. 13차는 점 14개를 모두 지나야 하니 잔차 제곱합이 0이어야 하는데 0.394가 나왔고, 11차의 0.861은 10차보다 큰 값이라 앞 절의 "차수를 올리면 커지지 않는다"에도 어긋납니다. 계산이 망가졌다는 증거입니다(망가진 숫자는 컴퓨터마다 조금 다를 수 있습니다).

지키는 법은 세 가지입니다.

정규방정식을 직접 풀지 않는다. np.linalg.lstsq(SVD)나 QR 분해는 $A$ 를 바로 다뤄 조건수를 제곱하지 않습니다. 실습대의 9차 맞춤도 브라우저 안에서 QR로 합니다.
$x$ 를 가운데로 옮기고 줄인다. $[-1, 1]$ 로 옮기면 조건수가 크게 줄어듭니다. Polynomial.fit이 이것을 자동으로 합니다.
차수를 낮게 둔다. 9-5에서 봤듯이 높은 차수는 계산 전에 모델로서도 나쁩니다.

책의 기호 읽는 법

책에 적힌 모양	읽는 법	뜻	파이썬에서는
$\hat{y}_i$	와이 햇 아이	모델이 예측한 값	`np.polyval(c, t)`
$r_i$ , $e_i$	잔차	관측 − 예측	`y - yhat`
$\bar{y}$	와이 바	평균	`y.mean()`
$S_r$ , SSE, RSS	잔차 제곱합	$\sum r_i^2$	`np.sum(r**2)`
$S_t$ , SST	총제곱합	$\sum (y_i - \bar{y})^2$	`np.sum((y - y.mean())**2)`
$R^2$ , $r^2$	결정계수	$1 - S_r/S_t$	`1 - SSR / SST`
$s_{y/x}$ , $s_e$	표준 오차	$\sqrt{S_r/(n-m-1)}$	`np.sqrt(SSR / (n - d - 1))`
$A^{\mathsf T}A\,\mathbf{c} = A^{\mathsf T}\mathbf{y}$	정규방정식	최소제곱 계수를 주는 연립방정식	`np.linalg.lstsq(A, y)`
$\arg\min_{a,b} S$	아그민	$S$ 를 가장 작게 하는 $a, b$	(최소제곱 함수들의 답)

실무에서는

np.polyfit(x, y, deg, w=None)            # 다항식 계수 (높은 차수부터), 가중치 w
np.polynomial.Polynomial.fit(x, y, deg)  # 새 방식: 정의역을 [-1, 1]로 옮겨 안정적, p(x)로 바로 계산
np.linalg.lstsq(A, y)                    # 열을 직접 고른 일반 최소제곱 (SVD)

scipy.stats.linregress(x, y)             # 직선 하나: 기울기, 절편, r, 기울기의 표준 오차
scipy.optimize.curve_fit(f, x, y, p0)    # 비선형 모델을 원래 공간에서 바로 맞춤 (로그 선형화의 함정 없이)

예측 구간(예: "20%에 닿는 시각은 34.5 ± 1분")까지 내야 한다면 statsmodels의 회귀 결과가 신뢰 구간을 계산해 줍니다. 기계 학습의 scikit-learn의 LinearRegression도 안에서 같은 최소제곱을 풉니다. 이상치가 걱정되면 제곱 대신 덜 민감한 손실을 쓰는 로버스트 회귀(scipy.optimize.least_squares의 loss="huber" 등)를 씁니다.

자주 하는 실수 (2분)

계수 순서 헷갈리기. np.polyfit은 높은 차수부터([기울기, 절편]), Polynomial은 낮은 차수부터 계수를 둡니다. 섞어 쓰면 기울기와 절편이 뒤바뀝니다.

>>> np.round(np.polynomial.Polynomial.fit(t, y, 1).convert().coef, 4)
array([101.12  ,  -2.2092])

R²만 보고 모델 고르기. 차수를 올리면 R²는 늘 오릅니다. 잔차의 무늬, 표준 오차가 줄어드는 폭, 떼어 둔 시험 점의 오차를 함께 보세요.

기록 밖을 아무렇지 않게 예측하기. 높은 차수일수록 기록 끝에서 곡선이 튑니다. 외삽은 짧게, 그리고 새 기록이 오면 다시 맞추세요.

0이나 음수에 로그 씌우기. 지수 모델 선형화는 $y > 0$일 때만 됩니다. 센서가 0을 찍으면 np.log가 -inf를 내고 맞춤 전체가 망가집니다.

연습 문제 (3분)

6호기의 다음 기록이 들어왔습니다. 28분에 37.4%입니다. 점 15개로 직선과 2차식을 다시 맞추고, 각각 20%에 닿는 시각을 구해 보세요.
9-3에서 직선의 잔차는 양 끝이 −, 가운데가 +였습니다. 직선과 2차식의 잔차에서 부호가 바뀌는 횟수를 세어 견줘 보세요. 잡음만 남았다면 점 14개에서 6~7번쯤 바뀌는 것이 보통입니다.
프로펠러 하나의 회전수(rpm)와 추력(N)을 재었습니다: 3000 → 5.1, 4000 → 9.0, 5000 → 14.2, 6000 → 20.3, 7000 → 27.8. 거듭제곱 모델 추력 = a·rpm^b 를 로그-로그 선형화로 맞추고 b를 구해 보세요.

풀이

>>> t15 = np.append(t, 28)
>>> y15 = np.append(y, 37.4)
>>> b1, a1 = np.polyfit(t15, y15, 1)
>>> round((20 - a1) / b1, 1)
np.float64(36.4)
>>> c15 = np.polyfit(t15, y15, 2)
>>> np.round(np.roots(c15 - [0, 0, 20]), 1)
array([-189.4,   34.5])
>>> def sign_changes(r):
...     return int(np.sum(np.diff(np.sign(r)) != 0))
>>> sign_changes(y - np.polyval(np.polyfit(t, y, 1), t)), sign_changes(y - np.polyval(c_2, t))
(2, 5)
>>> rpm = np.array([3000, 4000, 5000, 6000, 7000])
>>> thrust = np.array([5.1, 9.0, 14.2, 20.3, 27.8])
>>> b_p, lna_p = np.polyfit(np.log(rpm), np.log(thrust), 1)
>>> round(b_p, 3), round(np.exp(lna_p) * 1e7, 3)
(np.float64(2.002), np.float64(5.574))

1번은 새 점 하나로 직선의 예측이 36.7분에서 36.4분으로 당겨지고, 2차식은 34.5분 그대로입니다. 2차식이 이미 휘어짐을 잡고 있었기 때문에 새 기록에도 흔들리지 않습니다. 직선은 새 점이 들어올 때마다 조금씩 앞당겨지는데, 이것도 직선이 모델로 모자라다는 신호입니다.

2번은 직선의 잔차 부호가 두 번만 바뀝니다(− 네 개, + 여덟 개, − 두 개). 큰 무늬 하나가 있다는 뜻입니다. 2차식은 다섯 번 바뀌어 잡음에 훨씬 가깝습니다. 부호가 너무 드물게 바뀌면 모델이 무언가를 놓치고 있다는 신호입니다.

3번은 기울기 b가 약 2.00입니다. 추력은 회전수의 제곱에 비례한다는 프로펠러의 경험 법칙과 맞습니다. a는 아주 작은 수(약 5.57 × 10⁻⁷)라 보기 쉽게 10⁷을 곱해 적었습니다.

오 박사의 계산 노트

잔차 = 관측 − 예측. 최소제곱은 잔차 제곱합을 가장 작게 하는 계수를 고른다(매끄럽고 공식이 나온다, 대신 이상치에 약하다).
직선: b = (nΣty − ΣtΣy)/(nΣt² − (Σt)²), a = ȳ − b·t̄. 행렬로는 정규방정식 AᵀA c = Aᵀy.
평가: R² = 1 − S_r/S_t, 표준 오차 s = √(S_r/(n − 계수 개수)). R²보다 잔차의 무늬를 먼저 본다.
다항식·선형화도 같은 최소제곱. 로그로 바꾸면 로그 공간의 오차를 줄이니 가중치(w=y)로 보정한다.
차수를 올리면 훈련 오차는 줄지만 시험 오차와 외삽은 나빠진다(과적합). 기록 밖 예측은 짧게, 새 기록마다 다시.
NumPy: np.polyfit / np.polyval, np.polynomial.Polynomial.fit, np.linalg.lstsq. 정규방정식을 직접 풀지 않는다(조건수 제곱).

6호기는 마을 한 곳을 건너뛰고 돌아왔고, 착륙할 때 배터리는 21%였습니다. 다음 날 아침에는 민재 팀장이 돌풍 사고 보고서를 들고 옵니다. GPS 위치는 5초마다 찍혔는데, 돌풍은 12.3초에 불었습니다. 그 순간 드론이 어디 있었는지 알아야 합니다. 이번에는 잡음을 무시하고 큰 흐름을 그리는 회귀가 아니라, 기록된 점을 정확히 지나는 곡선이 필요합니다. 다음 교시는 보간과 스플라인입니다.

← 이전

8교시 · 지도가 커지면: LU 분해, 반복법, 조건수

10교시 · 비어 있는 시각 채우기: 보간과 스플라인

새 교시·새 강의 소식 받기

이 강의에 새 교시가 더해지거나 이어지는 강의가 열리면 이메일로 알려 드려요.