多项式拟合(polyfit)及局部加权回归(Lowess)对二维数据基础规律和离群特征学习的分析对比

概述：
1.当计算序列中离群区间的效应系数时，左右两侧增加的非离群区间应该尽量长些，离群程度越强，增加的非离群区间应越长，
多项式和Lowess才可能不被拉起。由于无法准确控制拉起的程度，则统一不拉起是更好的。
2.当需要学习序列的基本规律时，多项式可选尽量大的次数（在拟合自变量点上不会振荡），Lowess则是较大的frac和it。
解释：
对于多项式：
1.（特别重要）当多项式拟合次数过高时，新的应变量会在除了原自变量点之外的区间发生振荡，所以在应用拟合得到的该多项式时，新的自变量不能落在这些区间，
但若落在任意原自变量点上，则永远不会发生振荡。即虽然次数越高多项式越会过拟合，但若新自变量落在原自变量点上，则新应变量只会越来越趋近原应变量取值，
但永远不会振荡。
当拟合次数较低时，在所有区间和自变量点上都不易发生振荡，函数的跟随性较差，规律性较强。当新自变量落在原自变量点上时，新因变量也会偏离原因变量较远。
由此可知，当新应变量的取值不超过原自变量的集合时，拟合次数可适当取高，因为不用担心振荡；若拟合次数较低，新应变量会偏离原应变量更远，从而较不利于反映基础规律。
2.当拟合次数较高时，多项式容易在拟合时的自变量区间两端或一端及更外侧发生振荡，当取较中部区间的应变量时，是可以不发生振荡的。
3.当拟合点数过多，次数较高时，容易出现做最小二乘时奇异值分解不收敛的情况，从而无法拟合出多项式。
当拟合点数不多，拟合次数可以超过拟合点数，即未知数个数超过方程个数，虽然无法通过求解线性方程组得到未知数，但可通过梯度下降等最优化方法得到未知数，
即多项式的系数。此时容易发生振荡，中部非振荡区间也更窄。
对于局部加权回归：
1.Lowess方法返回的是拟合值，而不是函数，所以返回值与自变量x一一对应，而不能得到更密集或更稀疏或自变量区间外的返回值。
2.frac越接近1，拟合滑动窗口越接近全长，回归值越平滑，但frac不能大于1；it数越大，重新计算权重的最近点数目越多，受最近离群点影响越小；
因为此时已经比较平滑，所以所需it数较小。
3.frac越接近0，拟合滑动窗口中最近样本数越少，回归值跟随性越强，规律性越差；因为it数不宜大于滑动窗口中的样本数，所以此时it数也应较小。
4.当序列存在一段一段的离群值时，it值可取为最长离群段的样本数+1，这样无论用哪种Lowess方法，都不会受离群值的影响。
实践：

from math import ceil
import numpy as np
from scipy import linalg
import statsmodels.api as sm
import matplotlib.pyplot as plt
import copy
plt.style.use('seaborn-white')# Defining the bell shaped kernel function - used for plotting later on
def kernel_function(xi, x0, tau=.005):return np.exp(- (xi - x0) ** 2 / (2 * tau))def lowess_bell_shape_kern(x, y, tau=.005):"""lowess_bell_shape_kern(x, y, tau = .005) -> yestLocally weighted regression: fits a nonparametric regression curve to a scatterplot.The arrays x and y contain an equal number of elements; each pair(x[i], y[i]) defines a data point in the scatterplot. The function returnsthe estimated (smooth) values of y.The kernel function is the bell shaped function with parameter tau. Larger tau will result in asmoother curve."""n = len(x)yest = np.zeros(n)# Initializing all weights from the bell shape kernel functionw = np.array([np.exp(- (x - x[i]) ** 2 / (2 * tau)) for i in range(n)])# Looping through all x-pointsfor i in range(n):weights = w[:, i]b = np.array([np.sum(weights * y), np.sum(weights * y * x)])A = np.array([[np.sum(weights), np.sum(weights * x)],[np.sum(weights * x), np.sum(weights * x * x)]])theta = linalg.solve(A, b)yest[i] = theta[0] + theta[1] * x[i]return yestdef lowess_ag(x, y, frac=2/3, it=3):"""lowess(x, y, frac=2./3., it=3) -> yestLowess smoother: Robust locally weighted regression.The lowess function fits a nonparametric regression curve to a scatterplot.The arrays x and y contain an equal number of elements; each pair(x[i], y[i]) defines a data point in the scatterplot. The function returnsthe estimated (smooth) values of y.The smoothing span is given by frac. A larger value for frac will result in asmoother curve. The number of robustifying iterations is given by it. Thefunction will run faster with a smaller number of iterations."""n = len(x)r = int(ceil(frac * n))h = [np.sort(np.abs(x - x[i]))[r] for i in range(n)]w = np.clip(np.abs((x[:, None] - x[None, :]) / h), 0.0, 1.0)w = (1 - w ** 3) ** 3yest = np.zeros(n)delta = np.ones(n)for iteration in range(it):for i in range(n):weights = delta * w[:, i]b = np.array([np.sum(weights * y), np.sum(weights * y * x)])A = np.array([[np.sum(weights), np.sum(weights * x)],[np.sum(weights * x), np.sum(weights * x * x)]])beta = linalg.solve(A, b)yest[i] = beta[0] + beta[1] * x[i]residuals = y - yests = np.median(np.abs(residuals))delta = np.clip(residuals / (6.0 * s), -1, 1)delta = (1 - delta ** 2) ** 2return yestlowess_sm = sm.nonparametric.lowess
tau = 0.005  # abs(tau)越大，lowess_bell_shape_kern的回归值越平滑；tau→0+，回归值跟随性越强，规律性越差。
samples_in_window = 24  # samples_in_window越接近样本全体数量，回归值越平滑；越接近1，回归值跟随性越强，规律性越差。
frac = 0.9
n = 31  # 多项式最高次数，未知数（待拟合变量）有n+1个# 测试第一种正弦数据集
scale = 5
x_plot = np.random.uniform(low=0, high=1*np.pi, size=200*scale)
x_plot.sort()
y_plot = np.sin(x_plot * 1.5 * np.pi)
y_noise_plot = y_plot + np.random.normal(scale=1/2, size=len(x_plot))
x = np.array([x_plot[5*i] for i in range(int(len(x_plot)/scale))])
y = np.sin(x * 1.5 * np.pi)
y_noise = y + np.random.normal(scale=1/2, size=len(x))# lowess方法返回的是拟合值，而不是函数，所以返回值与自变量x一一对应，而不能得到更密集或更稀疏的返回值
yest_bell_plot = lowess_bell_shape_kern(x_plot, y_noise_plot, tau=tau)
# frac越接近1，拟合滑动窗口越接近全长，回归值越平滑；it数越大，重新计算权重的最近点数目越多，受最近离群点影响越小；因为此时已经比较平滑，所以所需it数较小。
# frac越接近0，拟合滑动窗口中最近样本数越少，回归值跟随性越强，规律性越差；因为it数不宜大于滑动窗口中的样本数，所以此时it数也应较小。
yest_ag_plot = lowess_ag(x_plot, y_noise_plot, frac=samples_in_window/len(x), it=3)
yest_sm_plot = lowess_sm(y_noise_plot, x_plot, frac=samples_in_window/len(x), it=3, return_sorted=False)
# polyfit与poly1d返回的是函数p，所以对这个函数p赋任何x，都可以得到相应的因变量y
c = np.polyfit(x, y_noise, n)
p = np.poly1d(c)
y_poly_plot = p(x_plot)# plot to compare three kinds of Lowess method
plt.figure(figsize=(20, 8))
plt.plot(x, y, color='g', label='sin()')
plt.scatter(x, y_noise, facecolors='none', edgecolor='darkblue', label='sin()+np.random.normal()')
plt.plot(x, y_noise, color='darkblue')
plt.fill(x[:], kernel_function(x[:], max(x)/2, tau), color='lime', alpha=0.5, label='Bell shape kernel')
plt.plot(x_plot, yest_ag_plot, color='orange', label='Lowess: A. Gramfort')
plt.plot(x_plot, yest_bell_plot, color='red', label='Lowess: bell shape kernel')
plt.plot(x_plot, yest_sm_plot, color='magenta', label='Lowess: statsmodel')
plt.plot(x_plot, y_poly_plot, color='k', label='polynomial of degree {}'.format(n))
plt.ylim(min(y_noise) - abs(min(y_noise) / 2), max(y_noise) * 1.5)  # 使y坐标的上下限关于y的最大最小值对称，且不受多项式拟合函数的振荡值影响
plt.legend()
plt.title('Sine with noise: Lowess regression comparing with polynomial')# 测试第二种带离群值的数据集
x = [np.linspace(1, 14, 14), np.linspace(1, 38, 38)]
# x_plot使绘图时的x取值比多项式拟合时，自变量x的取值更密集，以观察到可能出现的振荡
x_plot = [np.linspace(1, int(max(x[0])), int(max(x[0]))*5), np.linspace(1, int(max(x[1])), int(max(x[1]))*5)]
y = [np.random.randn(len(x[0])), np.random.randn(len(x[1]))]
# 设置离群值
for i in range(len(x)):y[i][int(len(x[i])/2-2):int(len(x[i])/2+2)] = 5*(i+1)
n = [[3, 4], [3, 6]]
weights = [(max(y[i])*2-y[i])/sum(max(y[i])*2-y[i]) for i in range(len(x))]  # y[i]越大，weights越小；让离群点占较小权重，正常点占较大权重。c = copy.deepcopy(n)  # c与n具有相同的结构；一定要用copy.deepcopy，不能用浅拷贝copy，否则后面对c的赋值也会改变n
p = copy.deepcopy(n)  # p与n具有相同的结构
yest_sm, yest_ag = [], []
for i in range(len(x)):for j in range(len(n[i])):c[i][j] = np.polyfit(x[i], y[i], n[i][j], full=True, w=weights[i])p[i][j] = np.poly1d(c[i][j][0])print('序列长度{0}个点，拟合次数为{1}，即目标函数未知数个数为{2}，即目标函数对自变量偏导数的方程个数为{2}：'.format(len(x[i]), n[i][j], n[i][j]+1))print('拟合函数的系数（次数由高到低）：', '\n', c[i][j][0])print('拟合值与实际值的MAPE：', sum(abs((p[i][j](x[i])-y[i])/y[i]))/len(x[i]))print('目标函数对自变量偏导数方程组的系数矩阵的秩：', c[i][j][2])print('目标函数对自变量偏导数方程组的系数矩阵的奇异值：', '\n', c[i][j][3])print('拟合的相关条件数：', c[i][j][4], '\n')plt.figure(figsize=(10, 6))plt.scatter(x[i], y[i], facecolors='none', edgecolor='darkblue', label='original_data')plt.plot(x[i], y[i], color='darkblue')plt.plot(x_plot[i], p[i][0](x_plot[i]), '-', color='g', label='polynomial of degree {}'.format(n[i][0]))plt.plot(x_plot[i], p[i][1](x_plot[i]), '--', color='orange', label='polynomial of degree {}'.format(n[i][1]))yest_sm.append(lowess_sm(y[i], x[i], frac=frac, it=5, return_sorted=False))yest_ag.append(lowess_ag(x[i], y[i], frac=frac, it=5))plt.plot(x[i], yest_sm[i], color='k', label='Lowess: statsmodel')plt.plot(x[i], yest_ag[i], color='r', label='Lowess: A. Gramfort')plt.ylim(min(y[i])-abs(min(y[i])/2), max(y[i])*1.5)  # 使y坐标的上下限关于y的最大最小值对称，且不受多项式拟合函数的振荡值影响plt.legend()plt.title('series length: {0}, degree of polyfit: {1}, {2}'.format(len(x[i]), n[i][0], n[i][1]))

图形展示：
sin()+noise：polynomial及Lowess均可较好地学习

两侧较短的非离群区间：polynomial及Lowess均被拉起

两侧较长的非离群区间：polynomial及Lowess均可不被拉起

多项式拟合(polyfit)及局部加权回归(Lowess)对二维数据基础规律和离群特征学习的分析对比相关推荐

python123英文字符的鲁棒_Robust Locally Weighted Regression 鲁棒局部加权回归 -R实现
鲁棒局部加权回归 Ljt 作为一个初学者,水平有限,欢迎交流指正. 算法参考文献: (1) Robust Locally Weighted Regression and Smoothing Scatt ...
R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测
全文下载链接:http://tecdat.cn/?p=22632 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的异常点进行建模的方法. 相关视频我们将对一种叫做STL的算法进行研究,STL是 ...
局部加权回归（LOWESS）
文章目录核函数叶帕涅奇尼科夫(epanechnikov)核函数立方核应用核函数局部加权回归增强鲁棒性双平方函数鲁棒局部加权线性回归局部加权线性回归:local weighted re ...
R语言平滑算法LOESS局部加权回归、三次样条、变化点检测拟合电视节目《白宫风云》在线收视率
最近我们被客户要求撰写关于LOESS的研究报告,包括一些图形和统计输出. 此示例基于电视节目的在线收视率.我们将从抓取数据开始. # 加载软件包. packages <- c("gpl ...
局部加权线性回归(Local Weighted Linear Regression)+局部加权回归+局部线性回归
局部加权线性回归(Local Weighted Linear Regression)+局部加权回归+局部线性回归 locally weighted scatterplot smoothing,LOWE ...
机器学习笔记：局部加权回归 LOESS
0 前言对于预测问题,回归中最简单的线性回归,是以线性的方法拟合出数据的趋势. 但是对于有周期性,波动性的数据,并不能简单以线性的方式拟合,否则模型会偏差较大局部加权回归(lowess)能较好的处 ...
STL——以鲁棒局部加权回归作为平滑方法的时间序列分解方法
摘要 STL是一种把时间序列分解为趋势项(trend component).季节项(seasonal component)和余项(remainder component)的过滤过程. STL有一个简单 ...
局部加权回归(LWR) Matlab模板
将百度文库上一份局部加权回归的代码,将其改为模板以便复用. q2x,q2y为数据集,是n*1的矩阵: r是波长参数,就是对于距离的惩罚力度: q_x是要拟合的数据横坐标,是1*n的矩阵: 得到的q_y ...
机器学习笔记（一）-局部加权回归（Locally weighted regression）LWR
在网上通过看斯坦福大学的机器学习课程,觉得讲的非常好.同时,为了加强自己的记忆,决定将自己学到的东西和一些理解记录下来,希望有所收获.废话不多说,直接开始笔记: 局部加权回归(locally weig ...

多项式拟合(polyfit)及局部加权回归(Lowess)对二维数据基础规律和离群特征学习的分析对比

多项式拟合(polyfit)及局部加权回归(Lowess)对二维数据基础规律和离群特征学习的分析对比相关推荐

最新文章

热门文章