크레이지J의 탐구생활

+최소제곱법

샘플데이터들을 최대한 근사할 수 있는 하나의 직선이나 평면으로 모델링하여 다른 샘플들도 모델링한 시스템으로 예측할 수 있다.

실제 해(샘플데이터)와 모델이 예측한 결과와의 차이가 최소가 되게하는 방식.

샘플데이터 (x, Y)

y=ax+b 인 모델로 예측. 

에러값을  오차^2의 합으로 한다.  (Y-(ax+b))^2 으로 해야 맞는데, 제곱이므로 상관없다.

J(Cost) =  

위 Cost Function을  a, b로 각각 편미분하여 0이 되도록 한다. 최적화 문제.  (위 cost function이 convex 형태이기 때문에 가능하다. 아니면 local minima 문제가 있음.)

위 식을 정리하면

 : 1번식 

   : 2번식

2식을 n으로 나누면

   : 3번식

3번식에서 b를 구해 번식에 대입한다.

a로 묶어주면

양변을 n으로 나눠주면

a를 뺀 왼쪽부분을 잘 보면 (제곱의 평균 - 평균의 제곱)이다.  어디서 많이 본 형태이다. 이것은 x의 분산이다.

분산은 또 편차 제곱의 평균이다. 즉, 왼쪽 부분을  다르게 표현하자면 아래와 같다... (전개해서 계산하면 똑같다는 것을 알 수 있다.)

양변에 n을 곱하면 

import numpy as np
import matplotlib.pyplot as plt

# some data
x=np.random.rand(7)*10
print('x=',x)

# black box + noise
y=x*1.5 + np.random.rand(len(x))*3
print('y=', y)

# model : y = ax + b
mx = np.average(x)
my = np.average(y)

print('mx=', mx, ' my=', my)
a1=np.average((x-mx)*(y-my))
a2=np.average(np.power((x-mx), 2))
a=a1/a2
b=my - a * mx
print('a=',a, 'b=',b)

plt.figure()
plt.scatter(x,y, marker='o', c='blue')

# predict
newx=np.arange(0, 10, 0.2)
yy = a*newx + b
print('predict yy=', yy)
plt.scatter(newx,yy, marker='+', c='red')

plt.show()

[결과]

x= [ 5.91242894  8.31159962  5.84424276  1.6816099   3.96831144  5.44434753
  4.99976856]
y= [ 11.3111561   14.27345181   9.25338631   4.42757815   7.11410915
   8.43423682   7.51826014]
mx= 5.16604410833  my= 8.90459692524
a= 1.49495983563 b= 1.18156847417

임의의 파란점을 기반으로 최소제곱법으로 주어진 모델로 추정선을 만들어 그림.

+ Projection / 프로젝션 / 투영

a벡터에 b벡터를 프로젝션시키면, b에서 a로 수선의 발을 내린다.  

project b onto a 라고 표현. (b를 a에 수직으로 내림. 수직으로 내리면 최단거리가 됨.)

수선의 발을 p라 하면. p점은 a 선상에 있고 p벡터를 a벡터의 상수배로 볼 수 있다.

b벡터 끝에서 p벡터에 가는 벡터와 a 벡터는 수직이다. 따라서 

수직인 벡터의 내적은 0이다. 열벡터에서 내적은 transpose의 곱으로 표현.

   ,     

여기서 매트릭스 연산식으로 바꾸고, 파라미터를 b로 보고 시스템을 꺼내보면. 프로젝션한 포인트 p는 

    ,    

이렇게 된다. 프로젝션 시스템 P는

이 변환 시스템이 어떤 점을 a 벡터로 프로젝션시키는 시스템이다.

프로젝션의 의미?

project b onto a (b를 a로 프로젝션 시킨다.)

     여기에서 는  a,b의 내적이다.  b를 a에 내린 사영의 크기 * a벡터의 크기이고,      는  a벡터 / a 벡터크기의 제곱.  즉, a 방향의 단위벡터 / a벡터의 크기가 되고,  이 둘을 곱하면  a벡터 크기 부분이 없어지므로  a 방향이고, 크기가 b를 a에 내린 사영의 크기인 벡터가 된다. 즉, 프로젝션이 된다.

b벡터에서 a벡터의 성분을 얼만큼 갖고있는지를 나타냄. 내적은 닮은 정도를 나타내기도 한다. 해당 성분이 얼만큼 있는지 나타내므로 닮은 정도이기도 하다. ( 내적이 a b cos (theta) 라는 의미를 생각해 보면 알 수 있다. 닮은 정도를 cos 유사도로 표현하는 것도 비슷하다. 두 벡터의 사이 각도가 0에 가까울 수록 비슷하고 cos값이 최대가 된다. 반대 방향으로 갈수록 작아진다.) 

+ 벡터의 투영과 최소 제곱

벡터와 평면이 수직인가를 알려면? 평면을  구성하는 기저벡터를 구하여 해당 벡터와 모두 수직인지 확인하면 된다.

임의 평면 x+2y+3z=0 의 기저 벡터를 구하려면?

위 방정식을 만족하는 독립인 벡터 두 개를 구하면 된다.

[1 1 1 ], [3 0 1]

위 두 개 벡터로 위 평면위의 모든 점을 스팬할 수 있다.

A는 위 기저벡터를 컬럼벡터로 하는 매트릭스로 어떤 시스템으로 볼 수 있다.

Ax는 위 평면위의 임의의 점이다. 기저벡터들의 선형결합으로 구성된 평면위의 점을 가리키는 벡터.

만약  평면 밖에 있는 b벡터를 A 평면에 내부의 벡터로 가장 가깝게 표현하려면??? 프로젝션!

즉, ||Ax-b||^2이 최소값이 되도록 한다. 최소제곱법.

Ax-b 벡터는 A 평면과 수직이다. 여기서 Ax는 b를 A 평면위에 프로젝션한 벡터이다. --> 최단 거리 지점.

b를 A에 프로젝션한 결과 : 

+ 실전

그렇다면 여러 점들을 가장 최적의 라인(또는 평면) 하나로 표현하려면?

측정값: (x1,y1), (x2,y2), (x3, y3), ...

피팅라인: y=ax+b   : 구하고자 하는 값 (a,b)

y1=ax1+b

y2=ax2+b

..

위 식을 모두 만족하는 a,b는 보통 존재하지 않는다. 다만 거리가 최소가 되도록 최적화된 라인을 찾을 수는 있다.

A x = b ( 왼쪽의 b는 위의 a,b 게수의 b와는 다른 것. y1,y2,y3를 말함)

x를 구하는 문제의 방정식을 만든다.

||Ax-b||^2 이 최소가 되게 하는 해를 구한다.

최소 거리를 만드는 문제는 프로젝션이다.

스페이스 밖의 점 b, 시스템은 A,   Ax는 b를 A에 프로젝션한 A 내부의 최적화된 점(p, 프로젝션)이 된다.

여기서 x를 구하면 최적화된 방정식의 계수인 a,b가 나온다.

예) (x1,x2) -> y    이러한 함수가 y=ax1+bx2+c 이 모델로 만들어진다고 한다. 아래 점들에 최대한 근사한  평면으로 표현하라면 최적의 계수 a,b,c?

(0,0) -> 1

(2,2)->6

(1,4)->2

(4,6)->10

A=np.array([[0,0,1],[2,2,1],[1,4,1],[4,6,1]])

b=np.array([[1],[6],[2],[10]])

np.dot( np.dot( np.linalg.inv( np.dot(np.transpose(A), A)), np.transpose(A) ) , b)

Out[72]: 

array([[ 3. ],

       [-0.5],

       [ 1. ]])

+ python source  

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


A=np.array([[0,0,1],[2,2,1],[1,4,1],[4,6,1]])
b=np.array([[1],[6],[2],[10]])

x = np.dot( np.dot( np.linalg.inv( np.dot(np.transpose(A), A)), np.transpose(A) ) , b)

print('A=',A)
print('b=',b)
print('x=',x)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(A[:,0], A[:,1], b[:,0], marker='o', label='sample')

ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('Y')
y = x[0]*A[:,0] + x[1]*A[:,1]+ x[2]
print('y=', y)
ax.scatter(A[:,0], A[:,1], y, marker='+', label='fitting')

# predict plane
x1=np.arange(0, 7, 0.2)
x2=np.arange(0, 7, 0.2)
x1,x2 = np.meshgrid(x1,x2)
Y = x[0]*x1 + x[1]*x2+ x[2]
ax.plot_surface(x1, x2, Y, alpha=0.2)


ax.legend()

plt.show()

파란점이 샘플 포인트이고 샘플 포인트에 최적화된 평면을 찾아 그 평면으로 다시 샘플포인트의 값을 추정하여 +로 표시하였음.   (거의 일치)