多项式拟合

多项式拟合

多项式的一般形式：

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()