如何设置HuberRegressor使0.1%样本为异常值及泛化方法问询

一、sklearn HuberRegressor的epsilon调整方法 #

sklearn的HuberRegressor中，异常值判定逻辑为：残差绝对值超过epsilon * sigma时被标记为异常，其中sigma是用**中位数绝对偏差（MAD）**估计的残差尺度。由于你的响应变量不服从正态分布，无法通过正态分布分位数直接计算epsilon，需通过迭代方式调整，让异常值占比逼近0.1%：

具体步骤 #

1. 初始拟合：先用普通线性回归得到初始残差，估计初始尺度sigma。
2. 计算目标阈值：根据当前残差分布，找到能让0.1%样本残差绝对值超过的阈值threshold。
3. 推导epsilon：epsilon = threshold / sigma（sigma为当前残差的MAD除以1.4826，该系数用于将MAD转换为与正态分布标准差可比的尺度）。
4. 迭代优化：用新epsilon拟合模型，重复步骤2-3，直到异常值占比稳定在目标范围内。

代码实现 #

import numpy as np
from sklearn.linear_model import HuberRegressor, LinearRegression

def find_huber_epsilon(X, y, target_outlier_ratio=0.001, max_iter=10):
    # 初始拟合用普通线性回归
    lr = LinearRegression(fit_intercept=True)
    lr.fit(X, y)
    y_pred = lr.predict(X)
    res = y.flatten() - y_pred
    
    for _ in range(max_iter):
        # 计算MAD并转换为尺度sigma
        mad = np.median(np.abs(res - np.median(res)))
        sigma = mad / 1.4826
        
        # 找到对应目标比例的残差阈值
        sorted_abs_res = np.sort(np.abs(res))
        threshold_index = int(len(sorted_abs_res) * (1 - target_outlier_ratio))
        threshold = sorted_abs_res[threshold_index]
        
        # 计算epsilon并重新拟合
        epsilon = threshold / sigma
        huber = HuberRegressor(epsilon=epsilon, fit_intercept=True, alpha=0)
        huber.fit(X, y)
        y_pred = huber.predict(X)
        res = y.flatten() - y_pred
        
        # 检查当前异常值比例，达标则提前终止
        current_ratio = huber.outliers_.sum() / len(y)
        if np.isclose(current_ratio, target_outlier_ratio, atol=1e-4):
            break
    
    return epsilon, huber

# 测试示例
X = np.random.rand(1000, 3)  # 放大样本量，0.1%比例更易观测
y = np.random.rand(1000, 1)
y[500] = 100  # 手动添加一个异常值

epsilon, model = find_huber_epsilon(X, y, target_outlier_ratio=0.001)
print(f"调整后的epsilon: {epsilon:.4f}")
print(f"异常值比例: {model.outliers_.sum() / len(y):.4f}")

二、通用方法适配其他带Huber损失的模型 #

对于GLS+Huber、GBRT+Huber这类模型，核心是根据目标异常值比例调整Huber损失中L1/L2损失的切换阈值，具体实现分两类：

1. 线性模型（如GLS+Huber） #

• 逻辑与sklearn HuberRegressor一致：
  • 初始化模型拟合得到残差，用MAD估计残差尺度。
  • 根据目标比例确定残差阈值，将其作为Huber损失的切换点。
  • 迭代拟合模型，更新残差和尺度，直到异常值比例符合要求。
• 可自定义Huber损失函数，结合scipy.optimize实现GLS+Huber的拟合，核心是在损失函数中加入残差阈值判断，切换L1/L2计算逻辑。

2. 树集成模型（如GBRT+Huber） #

多数树集成框架的Huber损失支持直接指定异常值比例：

• sklearn GradientBoostingRegressor：使用loss='huber'时，alpha参数对应“使用L2损失的样本比例”，设置alpha=0.999即可让0.1%样本作为异常值使用L1损失：

  from sklearn.ensemble import GradientBoostingRegressor
  gbrt = GradientBoostingRegressor(loss='huber', alpha=0.999, fit_intercept=True)
  gbrt.fit(X, y)
  

• LightGBM：objective='huber'支持huber_alpha参数，同样设置为0.999即可匹配0.1%的异常值比例。

内容的提问来源于stack exchange，提问作者Beitian Ma