使用Scipy进行线性回归的数值错误及标准误差疑问

我有两个数据数组：

• xdata：形状为(40,)，代表40年的平均海表温度
• ydata：形状为(40,721,1440)，维度依次为时间、纬度、经度，代表大气变量

我的目标是计算二者的回归斜率，并获取斜率误差以展示误差分布。尝试过程中遇到了这些问题：

1. 先通过计算协方差再开平方得到标准误差，为验证方法，调用scipy.stats.linregress时因数组形状不匹配触发ValueError: too many values to unpack (expected 4)
2. 改成循环遍历每个格点单独执行线性回归后，解决了函数调用问题，但得到的标准误差偏高，不确定是否正常，且期望误差呈现正态分布。

1. 搞定linregress的形状报错问题 #

scipy.stats.linregress默认只接受一维的x和y，直接传三维的ydata肯定报错。除了循环每个格点，用向量化操作能大幅提升效率，不用写嵌套循环：

• 先把xdata扩展维度，匹配ydata的时间轴：x = xdata[:, np.newaxis, np.newaxis]，形状变成(40,1,1)
• 要么用支持广播的新版scipy（≥1.7.0）直接调用linregress，要么手动实现向量化的线性回归计算，核心代码如下：

  import numpy as np
  
  # 计算均值
  x_mean = np.mean(xdata)
  y_mean = np.mean(ydata, axis=0)
  
  # 计算斜率的分子（协方差项）和分母（x的方差）
  numerator = np.sum((xdata - x_mean)[:, np.newaxis, np.newaxis] * (ydata - y_mean), axis=0)
  denominator = np.sum((xdata - x_mean)**2)
  
  # 回归斜率
  slope = numerator / denominator
  
  # 计算标准误差
  y_pred = slope * xdata[:, np.newaxis, np.newaxis] + (y_mean - slope * x_mean)
  residuals = ydata - y_pred
  n = len(xdata)
  se = np.sqrt(np.sum(residuals**2, axis=0) / (n-2)) / np.sqrt(np.sum((xdata - x_mean)**2))
  

2. 标准误差偏高的原因排查 #

循环得到的误差偏高不一定是错的，先排查这几个点：

• 数据本身特性：如果大气变量年际波动大，或者和海表温度的相关性弱，回归残差自然大，标准误差偏高是正常的
• 计算细节：检查循环里的残差自由度是不是用了n-2（线性回归拟合了斜率和截距两个参数，自由度要减2），用错成n-1的话误差会偏小，但如果是没减的话会偏大
• 缺失值影响：如果数据里有NaN，循环时没正确处理的话，有效样本量减少，也会导致标准误差变大

3. 验证误差的正态分布 #

线性回归里，斜率的抽样分布在满足高斯假设时是正态分布。可以这么验证：

• 把所有格点的标准误差拿出来画直方图，看形状是不是近似钟形
• 用Shapiro-Wilk或者Kolmogorov-Smirnov检验做统计验证，但注意大样本下，哪怕微小偏离也会让检验拒绝“正态分布”的原假设，所以主要看直方图的直观形态

手动计算和linregress的结果对比 #

可以随机挑几个格点，对比两种方法的结果，确认自己的计算没错：

from scipy.stats import linregress
import numpy as np

# 随机选一个格点
lat_idx = 100
lon_idx = 200
y_point = ydata[:, lat_idx, lon_idx]

# 用linregress计算
lr_slope, lr_intercept, lr_rvalue, lr_pvalue, lr_se = linregress(xdata, y_point)

# 手动计算
x_mean = np.mean(xdata)
y_mean = np.mean(y_point)
numerator = np.sum((xdata - x_mean)*(y_point - y_mean))
denominator = np.sum((xdata - x_mean)**2)
slope = numerator / denominator
residuals = y_point - (slope*xdata + (y_mean - slope*x_mean))
se = np.sqrt(np.sum(residuals**2)/(len(xdata)-2)) / np.sqrt(np.sum((xdata - x_mean)**2))

print(f"linregress斜率: {lr_slope}, 手动计算斜率: {slope}")
print(f"linregress标准误差: {lr_se}, 手动计算标准误差: {se}")

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