迭代拟合多项式曲线
Posted
技术标签:
【中文标题】迭代拟合多项式曲线【英文标题】：Iteratively fitting polynomial curve 【发布时间】：2019-09-05 00:30:03 【问题描述】：我想用以下方法迭代地拟合一条曲线到 python 中的数据：
-
拟合多项式曲线（或任何非线性方法）
丢弃值 > 曲线平均值的 2 个标准差
重复步骤 1 和 2，直到所有值都在曲线的置信区间内
我可以如下拟合多项式曲线：
vals = array([0.00441025, 0.0049001 , 0.01041189, 0.47368389, 0.34841961,
0.3487533 , 0.35067096, 0.31142986, 0.3268407 , 0.38099566,
0.3933048 , 0.3479948 , 0.02359819, 0.36329588, 0.42535543,
0.01308297, 0.53873956, 0.6511364 , 0.61865282, 0.64750302,
0.6630047 , 0.66744816, 0.71759617, 0.05965622, 0.71335208,
0.71992683, 0.61635697, 0.12985441, 0.73410642, 0.77318621,
0.75675988, 0.03003641, 0.77527201, 0.78673995, 0.05049178,
0.55139476, 0.02665514, 0.61664748, 0.81121749, 0.05521697,
0.63404375, 0.32649395, 0.36828268, 0.68981099, 0.02874863,
0.61574739])
x_values = np.linspace(0, 1, len(vals))
poly_degree = 3
coeffs = np.polyfit(x_values, vals, poly_degree)
poly_eqn = np.poly1d(coeffs)
y_hat = poly_eqn(x_values)
如何执行第 2 步和第 3 步？
【问题讨论】：
您只使用了 python 和 numpy 标签，您是否愿意使用其他 python 包，如 scipy、sklearn 等？ @JohnE，绝对，任何 Python 都可以 你能解释一下你的第二部分吗？特别是关于mean of curve。
【参考方案1】：
看起来你不会在这个过程中得到任何有价值的东西，有更好的技术来处理意外数据。谷歌搜索“异常值检测”将是一个好的开始。
话虽如此，下面是如何回答您的问题：
首先拉入库并获取一些数据：
import matplotlib.pyplot as plt
import numpy as np
Y = np.array([
0.00441025, 0.0049001 , 0.01041189, 0.47368389, 0.34841961,
0.3487533 , 0.35067096, 0.31142986, 0.3268407 , 0.38099566,
0.3933048 , 0.3479948 , 0.02359819, 0.36329588, 0.42535543,
0.01308297, 0.53873956, 0.6511364 , 0.61865282, 0.64750302,
0.6630047 , 0.66744816, 0.71759617, 0.05965622, 0.71335208,
0.71992683, 0.61635697, 0.12985441, 0.73410642, 0.77318621,
0.75675988, 0.03003641, 0.77527201, 0.78673995, 0.05049178,
0.55139476, 0.02665514, 0.61664748, 0.81121749, 0.05521697,
0.63404375, 0.32649395, 0.36828268, 0.68981099, 0.02874863,
0.61574739])
X = np.linspace(0, 1, len(Y))
接下来绘制数据的初始图：
plt.plot(X, Y, '.')
因为这可以让您了解我们正在处理的内容以及多项式是否适合 --- 简短的回答是，这种方法对于此类数据不会走得太远
此时我们应该停下来，但要回答这个问题，我将继续，主要遵循您的多项式拟合代码：
poly_degree = 5
sd_cutoff = 1 # 2 keeps everything
coeffs = np.polyfit(X, Y, poly_degree)
poly_eqn = np.poly1d(coeffs)
Y_hat = poly_eqn(X)
delta = Y - Y_hat
sd_p = np.std(delta)
ok = abs(delta) < sd_p * sd_cutoff
希望这是有道理的，我使用更高次的多项式，并且只在 1SD 处截断，否则不会丢弃任何东西。 ok 数组包含 True 值在 sd_cutoff 标准差内的那些点
为了检查这一点，我会再做一个情节。类似：
plt.scatter(X, Y, color=np.where(ok, 'k', 'r'))
plt.fill_between(
X,
Y_hat - sd_p * sd_cutoff,
Y_hat + sd_p * sd_cutoff,
color='#00000020')
plt.plot(X, Y_hat)
这给了我：
所以黑点是要保留的点（即X[ok] 给我这些点，np.where(ok) 给你指示）。
您可以使用这些参数，但您可能想要一个尾部较粗的分布（例如学生的 T 分布），但正如我上面所说，我的建议是使用 Google 进行异常值检测
【讨论】：
我看不出你是如何通过交替拟合和截面来实现第 3 点的...... @JirkaB。使用 2SD 的“截止”实际上不会删除任何内容，因此我的“2 保留所有内容”评论并没有实际实现循环。我刚刚注意到 OP 的评论，即“平滑方法可以更强大”，您的 RANSAC 建议将是一个不错的选择。我一直在使用贝叶斯高斯过程 MCMC 方法作为另一种选择，但它只有几百行代码，不确定是否值得发布。也许作为另一个答案，因为它可能对其他人很有趣 你甚至可以在 RANSAC 中集成贝叶斯高斯过程回归 您当然可以在 RANSAC 中放置一个 GP。称之为贝叶斯，我想整合所有参数的不确定性，scikit 实现似乎只是优化了超参数并忽略了任何不确定性。如果您只想要一个单点估计，那么这当然没问题，我只是更愿意对不确定性更诚实 这听起来很有希望......另一种方式可能是 ProSAC - researchgate.net/publication/4156175【参考方案2】：由于消除点与预期解决方案相距太远，您可能正在寻找RANSAC（随机抽样共识），即fitting a curve（或任何其他函数）到特定范围内的数据，例如您的情况 2 *性病。
您可以使用 scikit-learn RANSAC 估算器，它与包含的回归器（例如 LinearRegression）很好地对齐。对于您的多项式情况，您需要定义自己的回归类：
from sklearn.metrics import mean_squared_error
class PolynomialRegression(object):
def __init__(self, degree=3, coeffs=None):
self.degree = degree
self.coeffs = coeffs
def fit(self, X, y):
self.coeffs = np.polyfit(X.ravel(), y, self.degree)
def get_params(self, deep=False):
return 'coeffs': self.coeffs
def set_params(self, coeffs=None, random_state=None):
self.coeffs = coeffs
def predict(self, X):
poly_eqn = np.poly1d(self.coeffs)
y_hat = poly_eqn(X.ravel())
return y_hat
def score(self, X, y):
return mean_squared_error(y, self.predict(X))
然后你可以使用 RANSAC
from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(PolynomialRegression(degree=poly_degree),
residual_threshold=2 * np.std(y_vals),
random_state=0)
ransac.fit(np.expand_dims(x_vals, axis=1), y_vals)
inlier_mask = ransac.inlier_mask_
注意，X 变量被转换为 2d 数组，因为 sklearn RANSAC 实现需要它，并且在我们的自定义类中变平，因为 numpy polyfit 函数适用于 1d 数组。
y_hat = ransac.predict(np.expand_dims(x_vals, axis=1))
plt.plot(x_vals, y_vals, 'bx', label='input samples')
plt.plot(x_vals[inlier_mask], y_vals[inlier_mask], 'go', label='inliers (2*STD)')
plt.plot(x_vals, y_hat, 'r-', label='estimated curve')
此外，使用多项式阶数和剩余距离，我得到以下结果，度数=4，范围 1*STD
另一种选择是使用高阶回归器，例如Gaussian process
from sklearn.gaussian_process import GaussianProcessRegressor
ransac = RANSACRegressor(GaussianProcessRegressor(),
residual_threshold=np.std(y_vals))
说到对DataFrame的泛化，你只需要设置除一列之外的所有列都是特征，剩下的列是输出，如下所示：
import pandas as pd
df = pd.DataFrame(np.array([x_vals, y_vals]).T)
ransac.fit(df[df.columns[:-1]], df[df.columns[-1]])
y_hat = ransac.predict(df[df.columns[:-1]])
【讨论】：
感谢@Jirka B.！你知道是否可以完全忽略负异常值吗？ 负异常值是什么意思？一般来说，离群点是离预期太远的点，距离总是正的......【参考方案3】：需要三个函数来解决这个问题。首先需要一个线拟合函数来将一条线拟合到一组点：
def fit_line(x_values, vals, poly_degree):
coeffs = np.polyfit(x_values, vals, poly_degree)
poly_eqn = np.poly1d(coeffs)
y_hat = poly_eqn(x_values)
return poly_eqn, y_hat
我们需要知道从点到线的标准差。该函数计算标准差：
def compute_sd(x_values, vals, y_hat):
distances = []
for x,y, y1 in zip(x_values, vals, y_hat): distances.append(abs(y - y1))
return np.std(distances)
最后，我们需要比较一点到直线的距离。如果点到直线的距离大于标准差的两倍，则该点需要被剔除。
def compare_distances(x_values, vals):
new_vals, new_x_vals = [],[]
for x,y in zip(x_values, vals):
y1 = np.polyval(poly_eqn, x)
distance = abs(y - y1)
if distance < 2*sd:
plt.plot((x,x),(y,y1), c='g')
new_vals.append(y)
new_x_vals.append(x)
else:
plt.plot((x,x),(y,y1), c='r')
plt.scatter(x,y, c='r')
return new_vals, new_x_vals
如下图所示，此方法不适用于将线拟合到具有大量异常值的数据。由于距离拟合线太远，所有点最终都会被淘汰。
while len(vals)>0:
poly_eqn, y_hat = fit_line(x_values, vals, poly_degree)
plt.scatter(x_values, vals)
plt.plot(x_values, y_hat)
sd = compute_sd(x_values, vals, y_hat)
new_vals, new_x_vals = compare_distances(x_values, vals)
plt.show()
vals, x_values = np.array(new_vals), np.array(new_x_vals)
【讨论】：
以上是关于迭代拟合多项式曲线的主要内容，如果未能解决你的问题，请参考以下文章