如何在Python(Numpy)中应用1/σ²权重矩阵求解加权最小二乘

权重矩阵W的核心定义 #

加权最小二乘的权重矩阵W是对角矩阵，对角线上的每个元素对应单个样本的权重，即该样本噪声方差的倒数 1/σ²（σ为样本的噪声标准差）。噪声越大（σ越大）的样本，权重越小，拟合时对结果的影响也越小。

构造W的具体方法 #

根据你对噪声的了解，分两种常见场景处理：

• 场景1：已知每个样本的σ
  若你已经有每个样本的噪声标准差数组sigma，直接用np.diag()构造对角矩阵：

  W = np.diag(1 / sigma**2)
  

• 场景2：模拟/假设噪声分布
  比如你当前的代码中用了np.random.rand(50)生成噪声（均匀分布，范围[0,1]，方差为1/12），若假设所有样本σ相同，权重矩阵就是常数对角矩阵：

  # 均匀噪声的方差是1/12，所以1/σ²=12
  W = 12 * np.eye(50)
  

  若要模拟异方差（噪声随X变化），比如σ随X增大而增大，可先定义sigma数组再构造W：

  sigma = 0.1 + 0.5 * X  # 示例：σ随X线性增大
  W = np.diag(1 / sigma**2)
  

完整代码示例（含异方差模拟） #

import numpy as np
import matplotlib.pyplot as plt

# 生成带异方差噪声的样本数据
X = np.random.rand(50)
sigma = 0.1 + 0.5 * X  # 噪声标准差随X增大
Y = 2 + 3*X + np.random.normal(0, sigma, 50)  # 正态噪声贴合σ定义

plt.plot(X, Y, 'o')
plt.xlabel('X')
plt.ylabel('Y')

# 构造1/σ²权重矩阵
W = np.diag(1 / sigma**2)

# 构造含截距项的特征矩阵
X_b = np.c_[np.ones((50,1)), X]
# 计算加权最小二乘参数
beta = np.linalg.inv(X_b.T @ W @ X_b) @ X_b.T @ W @ Y

print(f"拟合截距: {beta[0]:.4f}, 拟合斜率: {beta[1]:.4f}")

# 绘制拟合直线
X_plot = np.linspace(0, 1, 100)
Y_plot = beta[0] + beta[1] * X_plot
plt.plot(X_plot, Y_plot, 'r-', label='加权拟合线')
plt.legend()
plt.show()

补充说明 #

• 用矩阵乘法@替代.dot()更符合Python的现代写法，两者功能一致
• 若样本量较大，直接求逆np.linalg.inv()可能效率低，可改用np.linalg.solve()求解线性方程组：

  beta = np.linalg.solve(X_b.T @ W @ X_b, X_b.T @ W @ Y)
  

  这种方法数值稳定性更好。

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