Python之数据分析（numpy中的多项式拟合）

向右看齐 2023-03-04 12:26 285阅读 0赞

**1、多项式拟合的概念**  
用一个无穷级数表示一个可微函数，任何可微的函数，总可以用一个N次多项式来近似，而比N次幂更高阶的部分可以作为无穷小量而被忽略不计。  
f(x) = p0x^n + p1x^n-1 + p2x^n-2 + … + pn  
y1 = f(x0)  
y2 = f(x1)  
…  
yn = f(xn)

**2、在numpy中求多项式拟合**

*  numpy.polyfit(自变量数组, 函数值数组, 最高次幂数n)  
    该方法将会返回\[p0, p1, …, pn\]多项式系数

**3、获取每一项的值**

*  numpy.polyval(\[p0, p1, …, pn\], 自变量数组)  
    返回多项式方程的根

**4、求微分**

*  numpy.polyder(\[p0, p1, …, pn\])  
    返回导函数系数数组

**5、案例分析**  
首先得到两曲线数据的价差：

diff_closing_price = apple_closing_prices - beer_closing_prices

然后拟合出原函数：

p = np.polyfit(days, diff_closing_price, 5)

之后取值得原函数的系数：

poly_closing_price = np.polyval(p, days)

若进行微分，则还可以求根、取实数根：

# 导函数的系数，即微分函数
    q = np.polyder(p)
    # 求方程的根，导函数方程的根，拐点
    roots = np.roots(q)
    # 根是取实根，因此需要判断一下，数组中哪些是实数
    reals = roots[np.isreal(roots)].real  # .real表示只取实部
    # print(reals)
    # 取实数中的每一个值
    peeks = [[days[0], np.polyval(p, days[0])]]
    for real in reals:
        if days[0]<real and real<days[-1]:
            peeks.append([real, np.polyval(p, real)])  # 拐点，和拐点对应的原函数值

**6、多项式拟合案例源码**

import datetime as dt
    import numpy as np
    import matplotlib.pylab as mp
    import matplotlib.dates as md
    
    
    def dmy2ymd(dmy):
        dmy = str(dmy, encoding='utf-8')  # 转码dmy日期
        date = dt.datetime.strptime(dmy, '%d-%m-%Y').date()  # 获取时间对象
        ymd = date.strftime('%Y-%m-%d')
        return ymd
    
    dates, beer_closing_prices = np.loadtxt(
        '0=数据源/beer_price2.csv', delimiter=',',
        usecols=(0, 4), unpack=True,
        dtype=np.dtype('M8[D], f8'),
        converters={ 0: dmy2ymd}
    )
    __, apple_closing_prices = np.loadtxt(
        '0=数据源/apple_price.csv', delimiter=',',
        usecols=(0, 4), unpack=True,
        dtype=np.dtype('M8[D], f8'),
        converters={ 0: dmy2ymd}
    )
    
    # 交易日两类数据的价差
    diff_closing_price = apple_closing_prices - beer_closing_prices
    # 将时间转化成以天为单位
    days = dates.astype(int)
    # 原函数
    p = np.polyfit(days, diff_closing_price, 5)
    # 5表示最高次数为5（一元五次函数），总共有6个方程（0次方也算一个）
    # 次数为1的时候就可以当做线性组合表示趋势：p = np.polyfit(days, diff_closing_price, 1)
    # 求多项式曲线的值，原函数的系数
    poly_closing_price = np.polyval(p, days)
    
    
    # 导函数的系数，即微分函数
    q = np.polyder(p)
    # 求方程的根，导函数方程的根，拐点
    roots = np.roots(q)
    # 根是取实根，因此需要判断一下，数组中哪些是实数
    reals = roots[np.isreal(roots)].real  # .real表示只取实部
    # print(reals)
    # 取实数中的每一个值
    peeks = [[days[0], np.polyval(p, days[0])]]
    for real in reals:
        if days[0]<real and real<days[-1]:
            peeks.append([real, np.polyval(p, real)])  # 拐点，和拐点对应的原函数值
    # 保证有序
    peeks.append([days[-1], np.polyval(p, days[-1])])
    peeks.sort()
    # 将列表peeks变成数组，二维数组，第一列是横坐标，第二列是纵坐标
    peeks = np.array(peeks)
    
    
    # 曲线图基础设置
    mp.figure('Polynomial Fitting', facecolor='lightgray')
    mp.title('Polynomial Fitting', fontsize=20)
    mp.xlabel('Date', fontsize=14)
    mp.ylabel('Difference Price', fontsize=14)
    
    # 主刻度设置为以周一为起始的星期格式
    ax = mp.gca()  # 获取刻度线（坐标轴）
    ax.xaxis.set_major_locator(
        md.WeekdayLocator(byweekday=md.MO)
    )
    # 次刻度设置为以天为单位
    ax.xaxis.set_minor_locator(
        md.DayLocator()
    )
    # 主刻度的格式化
    ax.xaxis.set_major_formatter(
        md.DateFormatter('%d %b %Y')
    )
    
    mp.tick_params(labelsize=10)  # 字体
    mp.grid(linestyle=':')  # 网格线
    
    # 绘制曲线
    dates = dates.astype(md.datetime.datetime)  # 将日期标准化成numpy的日期
    # 绘制多项式曲线
    mp.plot(dates, poly_closing_price, c="limegreen", linewidth=3, label="Polynomial Fitting")
    # 将每天的差价以散点的形式绘制出来
    mp.scatter(dates, diff_closing_price, c='dodgerblue', alpha=0.5, s=60, label='Difference Price')
    
    # 绘制导函数曲线，用箭头表示趋势
    dates, prices = np.hsplit(peeks, 2)  # 将peeks水平分割成两列（h）
    dates = dates.astype(int).astype('M8[D]').astype(md.datetime.datetime)
    # 这里分割出来的日期不是以天为单位的，需要转换成天，再转换成numpy专用的日期类型
    # 然后需要转化成绘图用的matplotlib的日期类型
    for i in range(1, dates.size):
        mp.annotate('', xytext=(dates[i-1], prices[i-1]),
                    xy=(dates[i], prices[i]),
                    size=40,
                    arrowprops=dict(arrowstyle='fancy', color="orangered", alpha=0.25)
                    )  # 添加文本注释
    # 画出散点
    mp.scatter(dates, prices, marker='^', color="orangered", s=80, label='Peek', zorder=4)
    
    mp.legend()  # 显示图例
    mp.gcf().autofmt_xdate()  # 设置格式展示的自动化调整
    mp.show()  # 显示图像

**7、拟合效果**  
![11][]

[11]: /images/20230208/21133a0aff5e4c0189285a926844c992.png