[Linear Analysis] 최소제곱법에 의한 수학적 모형

최소제곱법에 의한 수학적 모형

관련된 내용

• 최소제곱해

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+cx²
• 3차 다항식: y=a+bx+cx²+dx³

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+b{x}_1=\widehat{y_1},a+b{x}_2=\widehat{y_2},\cdots ,a+b{x}_n=\widehat{y_n}$$

위 식에서 각 식의 결과들을 y가 아닌 $\hat{y}$로 나타낸 것은 실험 결과에 의한 자료가 이론적으로 생성한 모형의 결과와 차이가 발생할 수 있기 때문입니다. 위 식을 행렬의 선형시스템의 형식으로 정리하여 보면 식 1과 같습니다.

$$\begin{array}{}\text{(식 1)}& \left[\begin{array}{rr}1& x_1\\ 1& x_2\\ ⋮& ⋮\\ 1& x_n\end{array}\right]\left[\begin{array}{r}a\\ b\end{array}\right]=\left[\begin{array}{r}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]\to Mv=y\end{array}$$

만약 위의 식을 $Mv=\hat{y}$로 나타냈다면 이 식은 자명한 해를 가지며 선형독립이 될 것입니다. 이것은 M을 구성하는 각 열벡터들은 $\hat{y}$에 대한 직교기저가 됩니다. 그러나 계산의 대상은 실제 값 y이므로 M의 열벡터들은 직교기저가 아닙니다. 즉, y와 $\hat{y}$ 사이에 오차가 존재하므로 정확한 해 대신 근사해를 계산해야 합니다. 이 계산은 식 2에 의해 이루어 집니다. (최소제곱해참조)

$$\begin{array}{rl}\text{(식 2)}& M^{T}Mv& =M^{T}yv& =(M^{T}M{)}^{-1}{M}^{T}y\end{array}$$

식 1과 2에 의해 생성된 식을 회귀식 또는 회귀모형(regression model)이라고 합니다. 실제값 y와 위 회귀모형에 의한 결과와 차이 즉, 오차를 잔차(residuals)라고 합니다(식 3).

$$\begin{array}{rl}\text{(식 3)}& d_1& =|y_1-(a+b{x}_1)|d_2& =|y_2-(a+b{x}_2)|\\ ⋮\\ d_n=|y_n-(a+b{x}_n)|\end{array}$$

식 3에서 d=잔차, 즉 회귀선의 값$\hat{y}$와 실제값 y 사이의 거리를 나타냅니다.

식 2가 성립하기 위해서는 $M^{T}M$가 역행렬이 존재하는 가역행렬이어야 합니다.

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^{T}M$은 가역행렬이므로 자명한 해를 가집니다. 즉, 모형의 계수가 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{array}{rl}\text{(식 4)}& & (x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)& y=a_0+a_{1}x+a_2{x}^2+\cdots +a_m{x}^m\end{array}$$

위 모형을 n개의 데이터에 적용하여 행렬형태로 정리하면 식 5와 같습니다.

$$\begin{array}{rl}\text{(식 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\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& x_n& x_n^2& \cdots & x_n^m\end{array}\right]\left[\begin{array}{r}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_m\end{array}\right]=\left[\begin{array}{r}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_n\end{array}\right]\to Mv=y& v=(M^{T}M{)}^{-1}{M}^{T}y\end{array}$$

예 2)

어떤 물체가 수직으로 하강할 경우 뉴턴의 제2 법칙을 따른다고 합니다.

$$\begin{array}{rl}s& =s_0+v_{0}t+\frac{1}{2}g{t}^2\\ s,s_0:& \text{하강거리, 시간(t)가 0일 때의 위치}\\ v_0,g:& \text{t=0일 때의 속도(v),중력가속도}\end{array}$$

실제 데이터가 다음과 같을 때 수학적 모형?

	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$$