手动计算回归beta估计量置信区间与内置函数结果不符排查

以下是可能导致手动计算结果与est.conf_int(alpha=alpha)偏差的关键问题，按优先级排序：

1. t值的分位数参数错误 #

这是最常见的错误来源：

• 代码中使用t_value = special.stdtrit(df, 1 - alpha)，但双侧置信区间需要取1 - alpha/2的分位数，而非1 - alpha。
• statsmodels的conf_int默认计算双侧置信区间，对应分位数为1 - alpha/2。例如alpha=0.05时，应取97.5%分位数，而非95%分位数。
• 修正方式：

  t_value = special.stdtrit(df, 1 - alpha/2)
  

2. X_estimator是否包含截距项 #

• 若X_estimator未添加全1的截距列，np.linalg.inv(np.transpose(X_estimator) @ X_estimator)的结果会与statsmodels内部的协方差矩阵（包含截距）不一致。
• 验证方法：检查X_estimator.shape[1]是否等于len(est.params)（即参数个数，包含截距）。若不等，需给X_estimator添加全1列：

  X_estimator = np.hstack([np.ones((X_estimator.shape[0], 1)), X_estimator])
  

3. 系数使用了四舍五入后的版本 #

• 代码中使用beta_estimator_4dec（四舍五入到4位小数的系数）进行计算，精度损失会导致置信区间偏差。
• 应直接使用statsmodels输出的原始系数：

  beta_estimator = est.params.values
  

4. 矩阵求逆的数值稳定性问题 #

• 直接使用np.linalg.inv对X.T @ X求逆可能存在数值误差，statsmodels内部使用更稳定的QR分解计算协方差矩阵。
• 替换为伪逆提升稳定性：

  V = np.linalg.pinv(np.transpose(X_estimator) @ X_estimator)
  

5. 自由度计算验证 #

• 手动计算的df = n - (k + 1)需与statsmodels的残差自由度一致：

  assert df == est.df_resid, "自由度计算不一致"
  

• 若不一致，检查k的定义：k应为不含截距的特征列数，确认X_train.columns是否包含截距项。

6. 标准误差的计算验证 #

• 对比手动计算的标准误与statsmodels输出的标准误：

  # statsmodels的标准误
  sm_std_errors = est.bse.values
  # 手动计算的标准误
  manual_std_errors = np.array(std_error_estimator_lst)
  # 查看差异
  print(np.abs(sm_std_errors - manual_std_errors))
  

• 若标准误本身存在偏差，优先排查X_estimator的定义和矩阵求逆步骤。

内容的提问来源于stack exchange，提问作者Maximilian Hauser