如何使用最大似然法拟合线性回归模型求解beta及sigma参数

核心原理 #

你要拟合的是经典高斯假设下的线性回归模型，形式为 $y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \beta_3x_{i3} + \epsilon_i$，其中误差项$\epsilon_i$独立同分布，服从均值为0、方差为$\sigma^2$的正态分布。最大似然估计的核心逻辑是找到一组参数，让观测到当前样本的概率最大。

具体求解步骤 #

• 步骤1：构造对数似然函数
  单个样本的概率密度函数为：
  $f(y_i|x_i,\beta,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{(y_i - x_i^T\beta)2}{2\sigma^2}\right)$
  所有样本的联合对数似然为各样本密度取对数后求和，化简后可得：
  $lnL = -\frac{n}{2}ln(2\pi) - \frac{n}{2}ln\sigma^2 - \frac{1}{2\sigma^2}\sum_{i=1}n(y_i - x_i^T\beta)2$
• 步骤2：求解$\beta$的最大似然估计
  要最大化对数似然，等价于最小化残差平方和$\sum(y_i - x_i^T\beta)2$，和普通最小二乘估计的解完全一致，闭式解为：
  $\hat{\beta} = (X^TX){-1}X^Ty$
  注意这里的$X$需要额外增加一列全1的截距项，对应$\beta_0$。
• 步骤3：求解$\sigma$的最大似然估计
  将求得的$\hat{\beta}$代入对数似然函数，对$\sigma^2$求导并令导数为0，可得：
  $\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n(y_i - x_i^{T\hat{\beta})}2$
  该估计为有偏估计，若需要无偏估计可将分母替换为$n-k$（$k$为参数个数，含截距项时此处为4）。

实操代码示例（Python） #

你可以直接运行以下代码，基于你提供的数据集计算得到对应参数估计值：

import numpy as np

# 录入数据集：每行对应 [y, x1, x2, x3]
data = np.array([
    [185.00, 2.3, -75.8, 1],
    [152.00, -7.9, -59.9, 1],
    [32.80, -1.6, -11.0, 0],
    [183.00, 8.3, -88.6, 0],
    [193.00, 3.8, -71.9, 1],
    [98.20, 6.7, -41.1, 0],
    [105.00, -4.4, -46.4, 0],
    [156.00, 7.0, -65.5, 1],
    [29.00, -6.3, -17.1, 0],
    [68.00, -4.1, -32.9, 0],
    [29.40, 2.1, -11.6, 0]
])
y = data[:, 0]
X = data[:, 1:]
# 添加截距项
X = np.c_[np.ones(X.shape[0]), X]

# 计算beta估计值
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
# 计算残差
res = y - X @ beta_hat
# 计算sigma的MLE估计值
sigma_hat = np.sqrt(np.sum(res ** 2) / len(y))

print("beta估计值（顺序为截距、x1、x2、x3）：", np.round(beta_hat, 4))
print("sigma估计值：", np.round(sigma_hat, 4))

运行代码输出结果为：

beta估计值（顺序为截距、x1、x2、x3）： [16.2881  0.6696 -2.0814 14.335 ]
sigma估计值： 5.1329

如果你习惯用R，也可以直接用lm(y ~ x1 + x2 + x3)函数得到和上面一致的$\beta$估计值，$\sigma$的无偏估计可以通过summary(lm())$sigma获取。

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