• 多项式拟合


    多项式拟合

    多项式的一般形式:

    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 = 4x3 + 3x2 - 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()

  • 相关阅读:
    Roads in the Kingdom CodeForces
    Vasya and Shifts CodeForces
    SOS dp
    Singer House CodeForces
    Codeforces Round #419 (Div. 1) (ABCD)
    在踏踏实实的生活里让自己坚持去做梦
    Replace Delegation with Inheritance
    Replace Inheritance with Delegation
    Replace Parameter with Methods
    Preserve Whole Object
  • 原文地址:https://www.cnblogs.com/maplethefox/p/11468296.html
Copyright © 2020-2023  润新知