Scikit-Learn HuberRegressor计算标准误、p值等统计量时奇异矩阵报错如何解决

报错原因 #

LinAlgError: Singular matrix 报错由两个核心问题导致：

1. 特征矩阵存在完全多重共线性，或样本数量≤特征维度，导致X^T X矩阵奇异，无法求普通逆
2. 直接套用OLS的标准误计算逻辑不适用于Huber这类稳健回归模型，即使解决奇异矩阵问题，输出的统计指标也不具备统计意义

解决方案 #

前置处理 #

先完成特征校验：

• 计算特征相关系数矩阵，删除完全共线的冗余特征
• 用VIF检验删除方差膨胀因子大于10的高共线性特征
• 确保样本数量远大于特征维度

完整实现代码 #

使用稳健三明治估计计算符合Huber回归假设的统计指标，同时用广义逆pinv替代普通逆inv兼容弱共线性场景：

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import HuberRegressor

# 1. 拟合Huber回归
solution = HuberRegressor()
solution.fit(beta_fac2, B)
params = np.append(solution.intercept_, solution.coef_)
predictions = solution.predict(beta_fac2)

# 2. 构造带截距项的特征矩阵
newX = np.append(np.ones((len(beta_fac2), 1)), beta_fac2, axis=1)
n, k = newX.shape  # n为样本量，k为参数个数

# 3. 计算Huber回归的残差和对应权重
residuals = B - predictions
huber_epsilon = solution.epsilon
# Huber权重规则：残差绝对值小于epsilon权重为1，否则为epsilon/|残差|
weights = np.where(np.abs(residuals) <= huber_epsilon, 1, huber_epsilon / np.abs(residuals))
Omega = np.diag(weights)

# 4. 计算稳健三明治协方差矩阵，用广义逆避免奇异矩阵问题
xtx = newX.T @ newX
xtx_inv = np.linalg.pinv(xtx)
xtox = newX.T @ Omega @ newX
cov_matrix = xtx_inv @ xtox @ xtx_inv
sd_b = np.sqrt(np.diag(cov_matrix))

# 5. 计算t统计量、p值、95%置信区间
ts_b = params / sd_b
df = n - k  # 自由度
p_values = [2 * (1 - stats.t.cdf(np.abs(i), df)) for i in ts_b]
t_critical = stats.t.ppf(0.975, df)
ci_low = params - t_critical * sd_b
ci_high = params + t_critical * sd_b

# 6. 整理输出结果
result = pd.DataFrame({
    "Coefficients": np.round(params, 4),
    "Standard Errors": np.round(sd_b, 3),
    "t values": np.round(ts_b, 3),
    "p values": np.round(p_values, 3),
    "95% CI Lower": np.round(ci_low, 4),
    "95% CI Upper": np.round(ci_high, 4)
}, index=["intercept"] + [f"feature_{i}" for i in range(beta_fac2.shape[1])])

print(result)

注意事项 #

• 如果替换为广义逆后仍存在数值不稳定问题，建议先对特征做标准化处理
• Huber回归的p值和置信区间属于近似推断结果，解读时建议结合多次抽样的稳健性检验

内容的提问来源于stack exchange，提问作者Albert Warnecki