多項式擬合

多項式擬合

多項式的一般形式：

y=p_{0}x^n + p_{1}x^{n-1} + p_{2}x^{n-2} + p_{3}x^{n-3} +...+p_{n}

多項式擬合的目的是為了找到一組p0-pn，使得擬合方程盡可能的與實際樣本數據相符合。

假設擬合得到的多項式如下：

f(x)=p_{0}x^n + p_{1}x^{n-1} + p_{2}x^{n-2} + p_{3}x^{n-3} +...+p_{n}

則擬合函數與真實結果的差方如下:

loss = (y_1-f(x_1))^2 + (y_2-f(x_2))^2 + ... + (y_n-f(x_n))^2

那么多項式擬合的過程即為求取一組p0-pn,使得loss的值最小。

X = [x1, x2, ..., xn] - 自變量
Y = [y1, y2, ..., yn] - 實際函數值
Y'= [y1',y2',...,yn'] - 擬合函數值
P = [p0, p1, ..., pn] - 多項式函數中的系數

根據一組樣本，並給出最高次冪，求出擬合系數
np.polyfit(X, Y, 最高次冪)->P

 

根據擬合系數與自變量求出擬合值, 由此可得擬合曲線坐標樣本數據 [X, Y']
np.polyval(P, X)->Y'

多項式函數求導，根據擬合系數求出多項式函數導函數的系數
np.polyder(P)->Q 

已知多項式系數Q 求多項式函數的根（與x軸交點的橫坐標）
xs = np.roots(Q)

兩個多項式函數的差函數的系數（可以通過差函數的根求取兩個曲線的交點）
Q = np.polysub(P1, P2)

案例：求多項式 y = 4x³ + 3x² - 1000x + 1曲線拐點的坐標。

'''
1. 求出多項式的導函數
2. 求出導函數的根，若導函數的根為實數，則該點則為曲線拐點。
'''
import numpy as np
import matplotlib.pyplot as mp

P = [4, 3, -1000, 1]
x = np.linspace(-20, 20, 1000)
# y = 4*x**3 + 3*x**2  - 1000*x + 1
y = np.polyval(P, x)  # 把x帶入P函數  得到y

# 求導
# Q = np.polyder([4,3,-1000,1])
Q = np.polyder(P)
xs = np.roots(Q)
# ys =  4*xs**3 + 3*xs**2  - 1000*xs + 1
ys = np.polyval(P, xs)
mp.plot(x, y)
mp.scatter(xs, ys, s=50, marker='o', c='orangered')
mp.show()

案例：使用多項式函數擬合兩只股票bhp、vale的差價函數：

# 多項式擬合
import numpy as np
import matplotlib.pyplot as mp
import datetime as dt
import matplotlib.dates as md


def dmy2ymd(dmy):
  """
  把日月年轉年月日
  :param day:
  :return:
  """
  dmy = str(dmy, encoding='utf-8')
  t = dt.datetime.strptime(dmy, '%d-%m-%Y')
  s = t.date().strftime('%Y-%m-%d')
  return s


dates, bhp_closing_prices = \
  np.loadtxt('bhp.csv',
             delimiter=',',
             usecols=(1, 6),
             unpack=True,
             dtype='M8[D],f8',
             converters={1: dmy2ymd})  # 日月年轉年月日
vale_closing_prices = \
  np.loadtxt('vale.csv',
             delimiter=',',
             usecols=(6,),
             unpack=True)  # 因為日期一樣,所以此處不讀日期
# print(dates)
# 繪制收盤價的折現圖
mp.figure('APPL', facecolor='lightgray')
mp.title('APPL', fontsize=18)
mp.xlabel('Date', fontsize=14)
mp.ylabel('Price', fontsize=14)
mp.grid(linestyle=":")

# 設置刻度定位器
# 每周一一個主刻度,一天一個次刻度

ax = mp.gca()
ma_loc = md.WeekdayLocator(byweekday=md.MO)
ax.xaxis.set_major_locator(ma_loc)
ax.xaxis.set_major_formatter(md.DateFormatter('%Y-%m-%d'))
ax.xaxis.set_minor_locator(md.DayLocator())
# 修改dates的dtype為md.datetime.datetiem
dates = dates.astype(md.datetime.datetime)

# 計算差價
diff_prices = bhp_closing_prices - vale_closing_prices
mp.plot(dates, diff_prices, color='dodgerblue', label='Diff Prices')

# 多項式擬合
days = dates.astype('M8[D]').astype('i4')
P = np.polyfit(days, diff_prices, 4)
y = np.polyval(P, days)
mp.plot(dates, y, color='orangered', linewidth=2, label='Polyfit line')

mp.legend()
mp.gcf().autofmt_xdate()
mp.show()