python做多项式拟合并绘图

最新推荐文章于 2024-11-12 10:40:54 发布

Hello AI！

最新推荐文章于 2024-11-12 10:40:54 发布

阅读量4k

点赞数 2

文章标签：深度学习 python 数据分析可视化数学建模

本文链接：https://blog.csdn.net/weixin_45081640/article/details/117717232

版权

本文所用文件的百度云链接：

链接：https://pan.baidu.com/s/15-qbrbtRs4frup24Y1i5og 
提取码：pm2c

之前有说过线性拟合了，显而易见，线性拟合在实际应用中局限性很大，多数时候并不能很好的描述数据的变换形势，这个时候就要考虑到使用非线性的方式，多项式拟合就是非线性拟合的其中一种方式，是相对简单的一种非线性的方式。

多项式拟合

多项式的一般形式:
$y = p_0x^n+p_1x^{n-1}+p_2x^{n-2}+...+p_n$
多项式拟合的目的是为了找到一组p₀-p_n , 使得拟合方程尽可能的与实际样本数据相符合.

假设拟合得到的多项式函数如下:
$f(x) = p_0x^n+p_1x^{n-1}+p_2x^{n-2}+...+p_n$
拟合得到的多项式函数与真实结果的误差方如下:
$loss=(y_1-f(x_1))^2 + (y_2-f(x_2))^2 + ...(y_n-f(x_n))^2$
那么多项式拟合的过程即为求取一组p₀-p_n , 使得loss的值最小.
这里的loss很重要，就是后面讲机器学习的时候，会说到损失函数，很多人对于损失函数的理解并不深刻，其实损失函数就是拟合结果与真实结果的偏差，就像这里的loss是方差的形势表现的，loss实际上是人为选择的，在这里也可以选择其他的方式，只是很多前辈们通过经验得出方差是描述多项式误差的一个很有效的数据，就像后面的机器学习所用到的损失函数，有很多也是直接调用API的，其实并不是一定要调用API，只是API调用的损失函数，是很多前辈们不断尝试之后，得到的最好的损失函数，比自己设计的效果好，并且节省时间。
多项式操作相关的API:

# 多项式拟合
# x, y:    为一组样本数据
# 最高次幂: 想得到的多项式函数的最高次幂  4
# p: 返回拟合得到的多项式的系数数组
p = np.polyfit(x,  y, 最高次幂) 
# 把x带入多项式函数p, 求得每个x对应的函数值
_y = np.polyval(p, x)
# 多项式求导  返回导函数的系数数组Q
Q = np.polyder(p)
# 已知多项式系数p, 求多项式函数的根(y=0时x的值)
xs = np.roots(p)

# 求两个多项式函数的差函数 (P1, P2)
Q = np.polysub(P1, P2)

案例:求多项式 y=4x³ + 3x² -1000x + 1 曲线的驻点坐标.
先来看下这里设计了一个一元三次的多项式，看下多项式的函数曲线。如图：

"""
多项式函数
"""
import numpy as np
import matplotlib.pyplot as mp

x = np.linspace(-20, 20, 1000)
y = 4*x**3 + 3*x**2 -1000*x + 1

# 求驻点坐标
P = np.array([4, 3, -1000, 1])
Q = np.polyder(P)
xs = np.roots(Q)
ys = np.polyval(P, xs)

mp.plot(x, y, color='dodgerblue')
mp.scatter(xs, ys, s=60, marker='s', c='red',
	zorder=3)
mp.show()

在这里插入图片描述

"""
多项式拟合
"""
import numpy as np
import matplotlib.pyplot as mp
import datetime as dt
import matplotlib.dates as md

# 当numpy解析文本时,将会把第一列中的每个字符串
# 都传给函数进行处理, 将处理完毕后的返回值
# 转成需要的M8[D]类型
def dmy2ymd(dmy):
	dmy = str(dmy, encoding='utf-8')
	# 把dmy转成日期对象
	d = dt.datetime.strptime(dmy, '%d-%m-%Y')
	t = d.date()
	s = t.strftime('%Y-%m-%d')
	return s

# 加载文件
dates, vale_closing_prices = np.loadtxt(
	'../da_data/vale.csv', delimiter=',', 
	usecols=(1,6), unpack=True, 
	dtype='M8[D], f8' , 
	converters={1:dmy2ymd})

bhp_closing_prices = np.loadtxt(
	'../da_data/bhp.csv', delimiter=',', 
	usecols=(6,), unpack=True)


# 绘制收盘价
mp.figure('DIFF', facecolor='lightgray')
mp.title('DIFF', fontsize=18)
mp.xlabel('Date', fontsize=14)
mp.ylabel('Price', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
# 设置主刻度定位器为每周一
ax = mp.gca()
ax.xaxis.set_major_locator(
	md.WeekdayLocator(byweekday=md.MO))
ax.xaxis.set_major_formatter(
	md.DateFormatter('%Y/%m/%d'))

# 把M8[D]转为matplotlib识别的date类型
dates = dates.astype(md.datetime.datetime)
diff = bhp_closing_prices - vale_closing_prices
mp.plot(dates, diff, 
	color='dodgerblue', linewidth=2,
	linestyle='--', label='diff line')

# 使用多项式拟合diff曲线
days = dates.astype('M8[D]').astype(int)
P = np.polyfit(days, diff, 4)
_y = np.polyval(P, days)
mp.plot(dates, _y, c='orangered', 
	linewidth=2, label='Polyfit Line')

mp.legend()
# 自动格式化x轴的日期输出
mp.gcf().autofmt_xdate()
mp.show()

重点看这几行代码


# 使用多项式拟合diff曲线
days = dates.astype('M8[D]').astype(int)
# 最高次幂: 想得到的多项式函数的最高次幂  4
P = np.polyfit(days, diff, 4)
# 把x带入多项式函数p, 求得每个x对应的函数值，也就是拟合的函数求出的y值
_y = np.polyval(P, days)
# 然后将这个拟合的多项式函数画出来
mp.plot(dates, _y, c='orangered', 
	linewidth=2, label='Polyfit Line')

在这里插入图片描述
如图所示：利用四次多项式拟合的曲线，拟合程度要远远高于线性函数，但是，这样也就又多了一个问题，那就是曲线过拟合。所以又要提到数据平滑处理，下面的平滑有用到卷积，其实原理很简单，就是做了一个卷积核，使突出的部分变的平滑，这样既可以防止过拟合，也可以平滑噪声的干扰。这个后面再细说，这里就简单说一下，主要是要明确拟合并不是和原始数据越契合就越好，有时候过拟合反而适得其反，所以既要排除噪声的干扰要防止过拟合不仅仅是做简单的数学曲线拟合，后面的机器学习中同样会有这样的问题。具体的处理方式，后面也会提到。

数据平滑

数据的平滑处理通常包含有降噪/拟合等操作.降噪的功能在于去除额外的影响因素. 拟合的目的在于数据的数学模型化, 可以通过更多的数学工具识别曲线特征.

案例: 绘制两只股票的收益率曲线.

收益率=(后一天收盘价-前一天收盘价) / 前一天收盘价

"""
demo06_ph.py  数据平滑处理
"""
import numpy as np
import matplotlib.pyplot as mp
import datetime as dt
import matplotlib.dates as md

# 当numpy解析文本时,将会把第一列中的每个字符串
# 都传给函数进行处理, 将处理完毕后的返回值
# 转成需要的M8[D]类型
def dmy2ymd(dmy):
	dmy = str(dmy, encoding='utf-8')
	# 把dmy转成日期对象
	d = dt.datetime.strptime(dmy, '%d-%m-%Y')
	t = d.date()
	s = t.strftime('%Y-%m-%d')
	return s

# 加载文件
dates, vale_closing_prices = np.loadtxt(
	'../da_data/vale.csv', delimiter=',', 
	usecols=(1,6), unpack=True, 
	dtype='M8[D], f8' , 
	converters={1:dmy2ymd})

bhp_closing_prices = np.loadtxt(
	'../da_data/bhp.csv', delimiter=',', 
	usecols=(6,), unpack=True)


# 绘制收盘价
mp.figure('Profit', facecolor='lightgray')
mp.title('Profit', fontsize=18)
mp.xlabel('Date', fontsize=14)
mp.ylabel('Price', fontsize=14)
mp.tick_params(labelsize=10)
mp.grid(linestyle=':')
# 设置主刻度定位器为每周一
ax = mp.gca()
ax.xaxis.set_major_locator(
	md.WeekdayLocator(byweekday=md.MO))
ax.xaxis.set_major_formatter(
	md.DateFormatter('%Y/%m/%d'))

# 把M8[D]转为matplotlib识别的date类型
dates = dates.astype(md.datetime.datetime)

# 整理两只股票的收益率
bhp_returns = np.diff(
	bhp_closing_prices) / bhp_closing_prices[:-1]
vale_returns = np.diff(
	vale_closing_prices) / vale_closing_prices[:-1]
dates = dates[:-1]

mp.plot(dates, bhp_returns, c='dodgerblue',
	label='bhp_returns',linewidth=1,
	alpha=0.5)
mp.plot(dates, vale_returns, c='orangered',
	label='vale_returns',linewidth=1,
	alpha=0.5)

# 卷积降噪
core = np.hanning(8)
core /= core.sum()
bhp_returns_convolved = \
	np.convolve(bhp_returns, core, 'valid')
vale_returns_convolved = \
	np.convolve(vale_returns, core, 'valid')
print(bhp_returns_convolved.size)

mp.plot(dates[7:], bhp_returns_convolved,
	c='dodgerblue', linewidth=2,
	label='bhp_convolved')
mp.plot(dates[7:], vale_returns_convolved,
	c='orangered', linewidth=2,
	label='vale_convolved')

mp.legend()
# 自动格式化x轴的日期输出
mp.gcf().autofmt_xdate()
mp.show()