logo标签列表

迭代拟合多项式曲线

pythonnumpyscikit-learnscipyscikits
9

我希望使用以下方法在Python中对数据进行迭代拟合曲线：

  1. 拟合多项式曲线（或任何非线性方法）
  2. 丢弃偏离曲线均值2个标准差以上的值
  3. 重复步骤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)

我该如何完成第二步和第三步？
你只使用了Python和NumPy标签，你是否愿意使用其他Python包，如SciPy、sklearn等？ - JohnE
@JohnE，完全没问题，任何Python都可以。 - user308827
你能再详细解释一下你的第二部分吗？特别是关于“曲线平均值”的部分。 - pouyan
3个回答
18
使用随机抽样一致性算法（RANSAC）来消除与期望解决方案相差太远的点，即对数据进行曲线拟合或其他函数拟合，并限制在特定范围内，例如您的情况下是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_

注意，由于sklearn RANSAC实现的要求，X变量被转换为二维数组，并在我们的自定义类中展开，因为numpy polyfit函数使用一维数组。

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')

visualisation of the poly-fitting

此外，通过调整多项式阶数和残差距离，我得到了以下结果，其中阶数为4，范围为1*STD。

enter image description here

另一个选择是使用高阶回归器，比如高斯过程。

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！你知道是否有可能完全忽略负的异常值吗？ - user308827
什么是负向异常值？通常来说，异常值是指距离期望太远的点，且这个距离总是正数... - Jirka
8

看起来按照那个步骤你不会得到任何有价值的东西，处理意外数据的方法有更好的技术。尝试搜索“异常检测”可能是一个好的开始。

话虽如此，这里是回答你问题的方法：

首先导入库并获取一些数据：

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)

这给了我：

data with poly and 1sd

所以黑点是需要保留的数据点（即X[ok]会将其返回，np.where(ok)会给出索引）。

您可以尝试调整参数，但可能需要一个具有更加粗尾巴的分布（例如学生T分布），但正如我上面所说，建议使用Google进行异常值检测。

我看不到你如何使用交替拟合和分段实现第三点... - Jirka
@JirkaB。使用“cutoff”为2SD并不会真正删除任何内容，因此我发表了“2保留所有内容”的评论，并没有实际实现循环。我刚刚注意到OP的评论，“平滑方法可以更加稳健”，您的RANSAC建议在这里是一个不错的选择。我一直在尝试使用贝叶斯高斯过程MCMC方法作为另一种选择，但它有几百行代码，不确定是否值得发布。也许作为另一个答案，因为它可能对其他人有趣。 - Sam Mason
你甚至可以在RANSAC中集成贝叶斯高斯过程回归。 - Jirka
你肯定可以在RANSAC中使用GP。如果要称其为贝叶斯方法，我希望能够对所有参数不确定性进行积分，而scikit实现似乎只是优化超参数并忽略任何不确定性。如果你只想要一个单点估计，那当然可以，但我更喜欢更诚实地表达不确定性。 - Sam Mason
听起来很有前途...另一种方法可能是ProSAC-https://www.researchgate.net/publication/4156175 - Jirka
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

如下图所示，当数据集中存在大量离群点时，该方法并不适合用于拟合一条直线。因为所有的点最终都会因为距离拟合直线太远而被排除掉。

elimination

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)
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接