Numpy의 polyfit과 poly1d의 사용법 - 최소제곱법과 polynomial class

제가 아주 예전에 공업수학 연재를 하면서 최소제곱법을 소개했던 적이 있습니다. 에러의 제곱의 합을 최소화하는 공업수학적 방법인데 아주 유용합니다. 그리고, 이를 이용한 Python의 Numpy 함수인 polyfit을 이용해서 최근 제가 집필한 책 파이썬으로 데이터 주무르기 1장에서 서울시 구별 CCTV의 수와 인구수와 관계를 직선으로 표현하려고 또 사용을 했죠. 초급자를 대상으로 해서, 머신러닝의 개념을 사용한 것은 아니었습니다. 그러다가 이 두 함수, polyfit과 poly1d의 사용예를 좀 더 보여드리고 싶다고 생각을 한거죠^^

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

t = np.arange(0, 10, 0.01)
y = 3*t + 5

plt.figure(figsize=(12,8))
plt.plot(t, y)
plt.show()

일단, 위의 코드로 시험 신호를 만들었습니다. 뭐 그냥 직선입니다.^^ 기울기가 3이고, 흔히 말하는 y 절편이 5입니다.^^

네. 저렇게 생긴거죠^^ 그리고, 저 직선을 찾는 재미(^^)를 주기 위해 살짝 노이즈를 한 두스푼(^^) 정도 넣습니다.

y_noise = y + np.random.randn(len(y))

plt.figure(figsize=(12,8))
plt.plot(t, y_noise)
plt.show()

그러면...

넵.. 저렇게 생긴거죠... 저게 만약 데이터고, 난 저 데이터에서 직선을 찾고 싶다는 거죠^^ 직선을 구성하는 것은 기울기와 절편입니다. 그 값을 알면 직선을 알게 되는 거죠...

fp1 = np.polyfit(t, y_noise, 1)
fp1

polyfit 명령으로 쉽게 찾을 수 있습니다. polyfit에서 세번쨰 인자는 찾고자 하는 함수의 차수입니다. 2를 넣어주면 2차식의 계수를 찾아달라고 하는 거죠.  아무튼 직선이니 저 결과는

저렇게 나옵니다. 애초 원 데이터의 기울기와 절편의 값과 유사합니다. 괜찮은거죠^^ 그런데.. 말이죠.. 저 계수를 가지고, 예측값을 찾을 수 없습니다. 함수(function)로 만들어 줘서 입력을 주고 결과를 얻을 수 있어야 하는거죠.

# poly1d : A one-dimensional polynomial class.
f1 = np.poly1d(fp1)
f1

위에 보이는 poly1d 함수를 사용해서 polynomial class를 만들어 주면 됩니다.

폴리노미얼 클래스는 수학정인 다항식으로 취급되면서 또 위 코드처럼 first class로 표현해서 코드에서 함수로 사용할 수도 있습니다. 여기서, polynomial 클래스 이야기를 조금 더 하면 말이죠~

위 두 일차식의 함을 표현해보면

np.poly1d([1, 1]) + np.poly1d([1, -1])

이렇게 됩니다. 이 결과는 당연히

이렇게 되죠... polynimial class가 이해되시나요? 하나 더???

그럼 곱하기 가야죠.. 우린 저 결과를 압니다만...

np.poly1d([1, 1]) * np.poly1d([1, -1])

해보면

이렇죠^^

이번엔.. 나눗셈 가보죠^^

np.poly1d([1, -2, 1]) / np.poly1d([1, -1])

이렇게 하면...

이렇게 되죠.. 바로.. 다항식 연산이 가능한 numpy의 클래스 중 하나입니다. 그리고, 오늘 이 poly1d 함수를 이용해서 우리도 결과를 얻어봐야죠.

plt.figure(figsize=(12,8))
plt.plot(t, y_noise, label='noise', color='y')
plt.plot(t, y, ls='dashed', lw=3, color='b', label='original')
plt.plot(t, f1(t), lw=2, color='r', label='polyfit')
plt.grid()
plt.legend()
plt.show()

위 코드를 보면, poly1d의 결과를 함수로 받은 f1이라는 아이에게 입력인자를 주는거죠. 그러면.. Python의 함수처럼 동작해서 결과를 줍니다.

짠... 어떤가요^^ 이번에는

y = np.square(t-1) + 1
y_noise = y + 10*np.random.randn(len(y))
fp1 = np.polyfit(t, y_noise, 2)
f1 = np.poly1d(fp1)

plt.figure(figsize=(12,8))
plt.plot(t, y_noise, label='noise', color='y')
plt.plot(t, y, ls='dashed', lw=3, color='b', label='original')
plt.plot(t, f1(t), lw=2, color='r', label='polyfit')
plt.grid()
plt.legend()
plt.show()

2차함수로 fit하고, 노이즈를 좀 많이 넣어보죠^^ 이제는 원 값과의 오차가 좀 생기겠죠^^

어때요? 괜찮죠...^^