최소제곱법에 의한 수학적 모형
관련된 내용
y=f(x)를 따르는 데이터에 대해 생각해 봅니다.
$$(x_1, y_2),\, (x_2, y_2), \, \cdot, \, (x_n, y_n)$$
위 데이터의 패턴을 찾는 것은 f(x)의 함수를 찾는 것과 같습니다. 이러한 함수를 수학적 모형이라 합니다. 이러한 함수의 몇가지 예를 나타내면 다음과 같습니다.
- 직선: y=ax+b
- 2차 다항식: y=a+bx+cx2
- 3차 다항식: y=a+bx+cx2+dx3
x1=symbols("x1")
f=2+3*x1+4*x1**2+5*x1**3
df=f.diff(x1)
ddf=df.diff(x1)
x=np.linspace(-10, 10, 100)
y1=[f.subs(x1, i) for i in x]
y2=[df.subs(x1, i) for i in x]
y3=[ddf.subs(x1, i) for i in x]
fig, ax=plt.subplots(figsize=(4,3))
ax.plot(x, y1, color="r", label=r"$y=dx^{3} + cx^{2} + b x + a$")
ax.plot(x, y2, color="b", label=r"$y=c x^{2} + b x + a$")
ax.plot(x, y3, color="g", label=r"$y=b x + a$")
ax.legend( bbox_to_anchor=(0.8,1))
ax.set_xlabel("x", loc="right")
ax.set_ylabel("y", loc="top")
ax.set_ylim(-200, 200)
ax.spines['left'].set_position(("data", 0))
ax.spines['bottom'].set_position(("data", 0))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()
실험적으로 획득한 다음의 데이터들에 대해 가장 간단한 모형인 직선을 고려해 봅니다.
$$(x_1, y_1), (x_2, y_2), \cdots , (x_n, y_n)$$
이 데이터들의 패턴으로 직선을 고려한다면 각 점들을 다음과 같이 나타낼 수 있습니다.
$$a+bx_1=\hat{y_1}, a+bx_2=\hat{y_2}, \cdots, a+bx_n=\hat{y_n}$$
위 식에서 각 식의 결과들을 y가 아닌 $\hat{y}$로 나타낸 것은 실험 결과에 의한 자료가 이론적으로 생성한 모형의 결과와 차이가 발생할 수 있기 때문입니다. 위 식을 행렬의 선형시스템의 형식으로 정리하여 보면 식 1과 같습니다.
$$\tag{식 1}\left[\begin{array}{rr}1&x_1\\1&x_2\\\vdots&\vdots\\1&x_n\end{array}\right]\left[\begin{array}{r}a\\b\end{array}\right]=\left[\begin{array}{r}y_1\\y_2\\\vdots\\y_n\end{array}\right] \rightarrow Mv=y$$
만약 위의 식을 $Mv=\hat{y}$로 나타냈다면 이 식은 자명한 해를 가지며 선형독립이 될 것입니다. 이것은 M을 구성하는 각 열벡터들은 $\hat{y}$에 대한 직교기저가 됩니다. 그러나 계산의 대상은 실제 값 y이므로 M의 열벡터들은 직교기저가 아닙니다. 즉, y와 $\hat{y}$ 사이에 오차가 존재하므로 정확한 해 대신 근사해를 계산해야 합니다. 이 계산은 식 2에 의해 이루어 집니다. (최소제곱해참조)
\begin{align} \tag{식 2}M^TMv&=M^Ty \\ v&=(M^TM)^{-1}M^Ty \end{align}
식 1과 2에 의해 생성된 식을 회귀식 또는 회귀모형(regression model)이라고 합니다. 실제값 y와 위 회귀모형에 의한 결과와 차이 즉, 오차를 잔차(residuals)라고 합니다(식 3).
\begin{align} \tag{식 3} d_1 &= |y_1 - (a+bx_1)|\\d_2 & = |y_2 - (a+bx_2)|\\ \qquad \vdots\\ d_n = |y_n - (a+bx_n)|\end{align}
식 3에서 d=잔차, 즉 회귀선의 값$\hat{y}$와 실제값 y 사이의 거리를 나타냅니다.
식 2가 성립하기 위해서는 $M^TM$가 역행렬이 존재하는 가역행렬이어야 합니다.
import numpy as np import numpy.linalg as LA from sympy import *
예 1)
(0, 1), (1, 3), (2, 4), (3, 4)이 4개의 점을 지나는 최소제곱직선(회귀직선)?
M=np.array([[1, 0], [1, 1],[1,2], [1,3]]) y=np.array([1, 3, 4, 4]).reshape(4, 1) MtM=np.dot(M.T, M) la.det(MtM).round(3)
20.0
$M^TM$은 가역행렬이므로 자명한 해를 가집니다. 즉, 모형의 계수가 1개씩만 존재합니다.
v=np.dot(np.dot(la.inv(MtM), M.T), y) print(v)
[[1.5] [1. ]]
위 결과를 정리하면 회귀선(최소제곱선)은 다음과 같습니다.
y=1.5 + x
예 2)
다음은 어떤 스프링이 추의 무게에 따라 늘어난 길이에 대한 자료입니다. 이 결과에 대한 1차 직선모형?
| weight(x) | 0 | 2 | 4 | 6 |
|---|---|---|---|---|
| length(y) | 6.1 | 7.6 | 8.7 | 10.4 |
M=np.array([[1, 0], [1, 2],[1,4], [1,6]]) y=np.array([6.1, 7.6, 8.7, 10.4]).reshape(4, 1) MtM=np.dot(M.T, M) la.det(MtM).round(3)
80.0
v=np.dot(np.dot(la.inv(MtM), M.T), y) print(v)
[[6.1] [0.7]]
결과적으로 length = 6.1+0.7weight 입니다.
직선에 대한 적용을 최소자승 다항식에도 동일하게 적용할 수 있습니다. n개의 데이터에 대한 수학적 모형을 다식 4와 같이 구축할 수 있습니다.
\begin{align} \tag{식 4}&(x_1, y_1), (x_2, y_2), \cdots , (x_n, y_n)\\ & y=a_0+a_1x+a_2x^2+ \cdots + a_mx^m\end{align}
위 모형을 n개의 데이터에 적용하여 행렬형태로 정리하면 식 5와 같습니다.
\begin{align} \tag{식 5}& \left[\begin{array}{rrrrr}1&x_1&x_1^2&\cdots&x_1^m\\1&x_2&x_2^2&\cdots&x_2^m\\\vdots&\vdots&\vdots&\vdots&\vdots\\1&x_n&x_n^2&\cdots&x_n^m\end{array}\right]\left[\begin{array}{r}a_0\\a_1\\a_2\\\vdots\\a_m\end{array}\right]=\left[\begin{array}{r}y_0\\y_1\\y_2\\\vdots\\y_n\end{array}\right] \rightarrow Mv=y\\ & v=(M^TM)^{-1}M^Ty\end{align}
예 2)
어떤 물체가 수직으로 하강할 경우 뉴턴의 제2 법칙을 따른다고 합니다.
\begin{align} s& =s_0+v_0t+\frac{1}{2}gt^2\\s, s_0:&\;\text{하강거리, 시간(t)가 0일 때의 위치}\\ v_0, g:&\;\text{t=0일 때의 속도(v),중력가속도}\end{align}
실제 데이터가 다음과 같을 때 수학적 모형?
| Time(sec, x) | distance(y) | |
|---|---|---|
| 0 | 1 | -0.18 |
| 1 | 2 | 0.31 |
| 2 | 3 | 1.03 |
| 3 | 4 | 2.48 |
| 4 | 5 | 3.73 |
s0=np.ones(5) t=np.array([1,2,3,4,5]) t2=t**2 M=np.c_[s0, t, t2] s=np.array([-0.18, 0.31, 1.03, 2.48, 3.73]).reshape(-1, 1) print(M)
[[ 1. 1. 1.] [ 1. 2. 4.] [ 1. 3. 9.] [ 1. 4. 16.] [ 1. 5. 25.]]
print(s)
[[-0.18] [ 0.31] [ 1.03] [ 2.48] [ 3.73]]
mtm=M.T@M la.det(mtm).round(3)
700.0
v=np.dot(np.dot(la.inv(mtm), M.T),y) print(v)
[[-0.398 ] [ 0.0347] [ 0.1607]]
위 결과에 의하면 다음과 같이 정리 할 수 있습니다.
$$s = -0.398+0.0347t+0.1607t^2$$
댓글
댓글 쓰기