如何基于最小二乘回归获取各系数对应的正态分布（Python实现）

首先明确核心逻辑：你提到的“均值为0的正态分布”是线性模型的观测误差项的分布，不是系数本身的分布。在经典最小二乘假设下，求解得到的每个系数a_i都服从正态分布，两个核心参数分别是均值（即lstsq返回的系数估计值本身）、方差（需要通过残差和输入矩阵额外计算，不是lstsq直接输出的结果）。

前置维度校验 #

先修正你描述里的维度匹配问题，否则矩阵运算会报错：

• 若输入矩阵M形状为(m, n)（m行n列，代表m组观测、n个特征），要满足M @ a = y的矩阵乘法规则，系数向量a的长度必须为n（形状为(n,)或(n,1)），输出向量y的长度必须为m（形状为(m,)或(m,1)），和你最初写的维度是反过来的，先调整数据维度再计算。

参数计算原理 #

• 系数a_i的正态分布均值：直接取scipy.linalg.lstsq(M, y)返回的解向量a的第i位元素即可
• 系数a_i的正态分布方差：是系数协方差矩阵的第i个对角线元素，协方差矩阵计算公式为：
  cov(a) = σ² * inverse(M.T @ M)
  其中σ²是观测误差的方差，用残差估计：σ² = 残差平方和 / (m - n)，m-n是模型自由度（样本量减特征数）。

可直接运行的Python代码 #

import numpy as np
from scipy.linalg import lstsq

# ----------------------
# 替换成你自己的真实数据
# 要求：M形状(m,n)，y形状(m,)
m = 120  # 样本量
n = 4    # 系数/特征个数
M = np.random.randn(m, n)  # 示例随机输入矩阵
true_a = np.array([1.5, -2, 0.3, 3])  # 模拟真实系数
y = M @ true_a + np.random.randn(m) * 0.5  # 模拟带噪声的观测y
# ----------------------

# 调用lstsq求解最小二乘
a, resids, rank, s = lstsq(M, y)

# 估计误差项方差σ²
if m <= n:
    raise ValueError("样本量小于等于特征数，为欠定系统，无法估计误差方差")
sigma_sq = np.sum(resids) / (m - n)

# 计算系数协方差矩阵
cov_a = sigma_sq * np.linalg.inv(M.T @ M)

# 提取每个系数的正态分布参数
for idx in range(n):
    a_mean = a[idx]
    a_std = np.sqrt(cov_a[idx, idx])
    a_var = cov_a[idx, idx]
    print(f"系数a[{idx}] 服从正态分布 N(μ={a_mean:.4f}, σ={a_std:.4f}, σ²={a_var:.4f})")

补充说明 #

• 如果不想手动计算协方差，可以直接用statsmodels库的OLS接口，拟合后会直接返回系数的标准差、置信区间等统计量，底层计算逻辑和上述代码完全一致
• 上述计算基于「观测误差独立同分布、服从均值为0的正态分布」的经典假设，如果你的数据存在异方差、序列相关等问题，需要调整协方差矩阵的计算方式（比如使用稳健标准误），属于进阶调整场景。

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