稳健线性估计方法比较

在数据分析和机器学习领域,经常需要对数据进行线性估计。然而,现实世界的数据往往包含噪声和异常值,这可能会影响估计的准确性。为了解决这个问题,研究者们开发了多种稳健的线性估计方法,如RANSAC、Theil-Sen和HuberRegressor。本文将介绍这些方法,并展示它们在拟合正弦函数时的表现。

首先,考虑一个简单的场景,即数据中只有模型误差,没有测量误差。在这种情况下,使用一个三次多项式来拟合正弦函数。可以看到,不同的稳健估计方法在这种情况下的表现。RANSAC在处理y方向上的强异常值时表现良好,而Theil-Sen则在处理X和y方向上的小异常值时表现更好。然而,当异常值较大时,Theil-Sen的性能会下降,甚至不如普通最小二乘法(OLS)。HuberRegressor虽然可以减轻异常值的影响,但它并不试图完全过滤掉异常值,因此在某些情况下,其性能可能不如RANSAC和Theil-Sen。

接下来,考虑数据中存在X方向的测量误差。在这种情况下,通过在数据中引入随机误差来模拟测量误差。可以看到,不同的稳健估计方法在处理X方向的测量误差时的表现。RANSAC和Theil-Sen在处理小异常值时表现良好,但在处理大异常值时,Theil-Sen的性能会下降。HuberRegressor在减轻异常值影响方面表现较好,但在完全过滤异常值方面不如RANSAC和Theil-Sen。

最后,考虑数据中存在y方向的测量误差。在这种情况下,通过在数据中引入随机误差来模拟测量误差。可以看到,不同的稳健估计方法在处理y方向的测量误差时的表现。RANSAC在处理y方向上的强异常值时表现良好,而Theil-Sen则在处理X和y方向上的小异常值时表现更好。然而,当异常值较大时,Theil-Sen的性能会下降,甚至不如普通最小二乘法(OLS)。HuberRegressor虽然可以减轻异常值的影响,但它并不试图完全过滤掉异常值,因此在某些情况下,其性能可能不如RANSAC和Theil-Sen。

为了评估不同稳健估计方法的性能,使用中位数绝对偏差(MAD)来衡量预测值与非污染数据之间的偏差。通过比较不同方法的MAD值,可以得出哪种方法在特定情况下表现最佳。

在实际应用中,可以根据数据的特点和异常值的类型来选择合适的稳健估计方法。例如,如果数据中存在大量的强异常值,可以选择RANSAC;如果数据中的异常值较小,可以选择Theil-Sen或HuberRegressor。此外,还可以通过调整这些方法的参数来进一步提高它们的性能。

代码示例

import numpy as np from matplotlib import pyplot as plt from sklearn.linear_model import HuberRegressor, LinearRegression, RANSACRegressor, TheilSenRegressor from sklearn.metrics import mean_squared_error from sklearn.pipeline import make_pipeline from sklearn.preprocessing import PolynomialFeatures np.random.seed(42) X = np.random.normal(size=400) y = np.sin(X) X = X[:, np.newaxis] X_test = np.random.normal(size=200) y_test = np.sin(X_test) X_test = X_test[:, np.newaxis] y_errors = y.copy() y_errors[::3] = 3 X_errors = X.copy() X_errors[::3] = 3 y_errors_large = y.copy() y_errors_large[::3] = 10 X_errors_large = X.copy() X_errors_large[::3] = 10 estimators = [ ("OLS", LinearRegression()), ("Theil-Sen", TheilSenRegressor(random_state=42)), ("RANSAC", RANSACRegressor(random_state=42)), ("HuberRegressor", HuberRegressor()), ] colors = { "OLS": "turquoise", "Theil-Sen": "gold", "RANSAC": "lightgreen", "HuberRegressor": "black", } linestyle = { "OLS": "-", "Theil-Sen": "-.", "RANSAC": "--", "HuberRegressor": "--" } lw = 3 x_plot = np.linspace(X.min(), X.max()) for title, this_X, this_y in [ ("Modeling Errors Only", X, y), ("Corrupt X, Small Deviants", X_errors, y), ("Corrupt y, Small Deviants", X, y_errors), ("Corrupt X, Large Deviants", X_errors_large, y), ("Corrupt y, Large Deviants", X, y_errors_large), ]: plt.figure(figsize=(5, 4)) plt.plot(this_X[:, 0], this_y, "b+") for name, estimator in estimators: model = make_pipeline(PolynomialFeatures(3), estimator) model.fit(this_X, this_y) mse = mean_squared_error(model.predict(X_test), y_test) y_plot = model.predict(x_plot[:, np.newaxis]) plt.plot(x_plot, y_plot, color=colors[name], linestyle=linestyle[name], linewidth=lw, label="%s: error = %.3f" % (name, mse)) legend_title = "Error of Mean\nAbsolute Deviation\nto Non-corrupt Data" legend = plt.legend(loc="upper right", frameon=False, title=legend_title, prop=dict(size="x-small")) plt.xlim(-4, 10.2) plt.ylim(-2, 10.2) plt.title(title) plt.show()
沪ICP备2024098111号-1
上海秋旦网络科技中心:上海市奉贤区金大公路8218号1幢 联系电话:17898875485