基于Scipy与Sklearn的回归统计计算正确性及误差估计问询

1. MSE与MAE（你代码中写的MUE应为MAE）的误差计算 #

你当前使用的np.sqrt(2 / len(values_array))公式不准确：

• MSE的标准误：MSE是模型残差方差的无偏估计，其标准误的正确解析公式为：

  n = len(values_array)
  df = n - 2  # 线性回归的自由度（斜率+截距）
  mse_se = mse * np.sqrt(2 / df)
  

  这是基于残差服从正态分布的假设，方差为2*σ⁴/(n-k)（k为模型参数数量），用样本MSE代替σ²得到的近似。
• MAE的标准误：没有简单的解析公式，你当前的公式完全不适用。推荐用bootstrap方法估计：

  def bootstrap_mae_se(x, y, slope, intercept, n_bootstraps=1000):
      mae_list = []
      for _ in range(n_bootstraps):
          indices = np.random.choice(range(len(x)), len(x), replace=True)
          boot_x = x[indices]
          boot_y = y[indices]
          pred_y = slope * boot_x + intercept
          mae_list.append(mean_absolute_error(boot_y, pred_y))
      return np.std(mae_list)
  

2. 相关系数的手动实现与标准误 #

Pearson系数 #

你的手动实现存在明显bug：y_mean = np.mean(values_array)应该改为y_mean = np.mean(ordered_values_array)，分子部分也错误地用了(values_array - x_mean) * (values_array - y_mean)，正确的分子是(values_array - x_mean) * (ordered_values_array - y_mean)。修正后：

x_mean = np.mean(values_array)
y_mean = np.mean(ordered_values_array)
pearson_numerator = np.sum((values_array - x_mean) * (ordered_values_array - y_mean))
pearson_denominator = np.sqrt(np.sum((values_array - x_mean)**2) * np.sum((ordered_values_array - y_mean)**2))
pearson_r_manual = pearson_numerator / pearson_denominator

标准误公式np.sqrt((1 - pearson_r_manual ** 2) / (len(values_array) - 2))是正确的（基于双变量正态假设）。

Kendall系数 #

你的手动实现未处理数据中的结（相同值），而scipy.stats.kendalltau会自动调整结的影响，因此当数据存在结时，手动计算的τ会与scipy结果不一致。
标准误公式np.sqrt((2 * (2 * len(values_array) + 5)) / (9 * len(values_array) * (len(values_array) - 1)))仅适用于无结且样本量较大的情况，有结时需要使用调整后的公式，或直接用scipy返回的p值对应的标准误（可通过z = tau / se反推）。

Spearman系数 #

你的手动实现同样未处理结，当数据存在结时，scipy的spearmanr会使用调整后的公式（比如修正秩的 Pearson 相关），而你的公式仅适用于无结场景。

3. Spearman系数的Bootstrap标准误 #

你的Bootstrap流程逻辑正确，但有两个优化点：

• n_bootstraps=100样本量太小，建议至少设置为1000次，提升结果稳定性；
• 解析方法补充：当样本量较大且无结时，可使用近似公式se = np.sqrt((1 - rho**2)/(n-2))（与Pearson的标准误公式一致，因为Spearman rho是秩的Pearson相关），但Bootstrap方法在有结、样本量小的场景下更可靠。

优化后的Bootstrap函数：

def bootstrap_spearman(x, y, n_bootstraps=1000):
    bootstrapped_rhos = []
    for _ in range(n_bootstraps):
        indices = np.random.choice(range(len(x)), len(x), replace=True)
        boot_x = x[indices]
        boot_y = y[indices]
        rho, _ = spearmanr(boot_x, boot_y)
        bootstrapped_rhos.append(rho)
    return np.std(bootstrapped_rhos)

4. 斜率与截距的置信区间 #

你的截距标准误计算错误，正确的截距标准误公式为：

x = ordered_values_array
x_mean = np.mean(x)
ssxx = np.sum((x - x_mean)**2)
stderr_intercept = std_err * np.sqrt(1/len(x) + (x_mean**2)/ssxx)

斜率的标准误stderr_slope = std_err是正确的（linregress返回的std_err就是斜率的标准误），因此斜率的置信区间计算正确，但截距的置信区间需要用修正后的stderr_intercept重新计算。

验证数据集与预期结果 #

生成一个可控的线性数据集用于测试：

# 生成验证数据
np.random.seed(42)
x = np.linspace(0, 10, 50)
y_true = 2 * x + 3
y = y_true + np.random.normal(0, 1.5, size=len(x))  # 加入正态噪声，标准差1.5

# 预期指标参考
# Pearson r: ~0.98（高度线性相关）
# Spearman rho: ~0.98（秩相关接近线性相关）
# Kendall tau: ~0.92（秩一致性高）
# MSE: ~2.25（噪声方差的近似）
# MAE: ~1.2（正态分布下MAE≈σ*sqrt(2/π)≈1.5*0.797≈1.196）
# 斜率: ~2，截距: ~3
# 斜率95%置信区间: 2±t(48)*std_err（t(48,0.975)≈2.01，std_err≈0.05，区间约(1.9,2.1)）

将你的代码中的values_array替换为y，ordered_values_array替换为x，即可验证各指标计算的正确性。

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