在数据分析和机器学习领域,经常需要对数据进行线性估计。然而,现实世界的数据往往包含噪声和异常值,这可能会影响估计的准确性。为了解决这个问题,研究者们开发了多种稳健的线性估计方法,如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()