迭代拟合多项式曲线

Question

我想用以下方法迭代拟合曲线到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)

如何执行步骤2和3？

Answer 1

由于消除点离预期的解决方案太远，您可能正在寻找RANSAC（随机抽样共识），它可以将曲线（或任何其他函数）拟合到一定范围内的数据，例如您使用2 * STD的情况。

您可以使用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变量将转换为2d数组；在我们的自定义类中，由于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

另一种选择是使用像高斯过程这样的高阶回归

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

Answer 2

在执行该过程后，您似乎一无所获，有更好的技术来处理意外数据。谷歌搜索“异常检测”将是一个好的开始。

这样说，这是如何回答您的问题：

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

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进行离群值检测是我的建议