机器算法验证 - 通过矩阵方程手动拟合的简单线性回归与 lm() 输出不匹配 - 吾爱随笔录

通过矩阵方程手动拟合的简单线性回归与 lm() 输出不匹配

机器算法验证 r 回归自习线性代数

2022-03-22 12:05:05

我正在尝试使用矩阵将线性模型拟合到我的数据集，即使我可以使用 OLS 并且在没有矩阵的情况下执行此操作作为一个简单的教程，让我自己更好地理解R矩阵表示法。

这是我想要拟合的模型：

Y = X β + ε

$\bf Y=X\boldsymbol\beta+\varepsilon$

其中 $\bf Y$ 是一个 $1\times n$ 矩阵， $\bf X$ 是一个 $n\times k$ 矩阵（其中 $k$ 是 $\beta$ 的数量，在这种情况下为 2）， $\boldsymbol \beta$ 是一个 $k\times 1$ 矩阵，最后我们的误差项是 $n\times 1$ 。我理解这部分。

当我简单地使用lm()命令来拟合我的数据时，我从命令中得到以下信息summary()：

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0503 -1.4390  0.4921  1.0589  3.9446 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9849     0.8219   3.632  0.00191 ** 
x             0.5612     0.1084   5.178 6.32e-05 ***

所以summary()告诉我 $\beta$ 矩阵是一个 $2\times 1$ 矩阵，第一个数字（即 $\beta_0$ ）为2.0949，第二个数字（ $\beta_1)$ 为0.1084。我的问题是这样的：

我们知道矩阵 $\beta$ 实际上是：

β = (X^{T} X)^{- 1} (X^{T} Y)

$\boldsymbol\beta=(\bf{X}^T\bf{X})^{-1}(\bf{X}^T\bf{Y})$

当我尝试使用 R using 简单地手动执行此计算时b=(t(x)*x)^-1*t(x)*y，我得到一个 $1\times 20$ 向量（其中 $20$ 当然是 $n$ ，观察次数）。为什么我没有得到我应该得到的 $2\times 1$ 矩阵？

3个回答

您在 R 代码中犯了两个错误b。

solve用于矩阵求逆。提高X到幂会反转的每个元素，这有时会很有用，但这不是我们想要的。 $-1$ X
R 使用运算符%*%进行矩阵乘法。否则，它会进行元素乘法，并要求您的数组符合 R 对向量的处理。再次，偶尔有用。

这一切都在 R 文档中进行了更详细的解释，或者在线 R 材料的几乎不可数的介绍，例如这个。

b<-solve(t(X)%*%X)%*%t(X)%*%y是正规方程的字面表示，而b<-solve(crossprod(X), crossprod(X,y))速度更快，更惯用，并且与的输出相匹配lm()。但是您绝对不想出于数字原因直接显式计算正规方程（并且在您丢失所有有效数字之前，R 不会警告您）。该lm方法默认使用 QR 分解，在数值上更稳定。如果您应用正确的“手摇”计算并且它仍然与 R 输出不同，这可能是因为的显式反转导致数值精度损失。 $X^TX$

不要忘记为截距 $1$

您需要在矩阵中包含一列x，对应于截距。该lm()函数会自动为您执行此操作，但您需要在使用正规方程计算答案时自行添加。

请参阅本网站上的上一个问题：使用正规方程计算多元线性回归中的系数。有示例代码，您应该可以直接复制和粘贴。

这是一本很棒的入门书，解释了为什么必须添加列：http: //www.amazon.com/Applied-Regression-Models-Edition-Student/dp/0073014664

您还需要遵循@user777 给出的建议

此 R 代码可用于通过矩阵回归计算我称为的给定数据集的 Y（y 值的向量，拟合值）和 Beta（系数的向量）insert.dataset。即使您向公式中添加其他数字变量，这也应该有效。

library(matlib) # enables function inv() to calculate a matrix's inverse

model <- lm(formula = y ~ x, data = insert.dataset) # save linear model
beta0 <- rep(1, nrow(model$model)) # column of 1s representing coefficient beta0 (intercept)
X <- as.matrix(cbind(beta0, model$model[,-1]), nrow=nrow(model$model)) # create X matrix, replacing column of outcomes with beta0

# Matrix equation to create Y (fitted values), using X and coefficients
Y <- X %*% model$coefficients
model$fitted.values # should be identical to Y

# Matrix equation to create Beta (coefficients), using X and Y
Beta <- inv(t(X) %*% X) %*% t(X) %*% model$fitted.values
model$coefficients # should be identical to Beta

应该注意的是，来自的输出lm已经有表示 Y(model$fitted.values)和 Beta的向量model$coefficients，但是需要进行一些小的修改才能获得 X（观察到的预测值的矩阵）。这个小的修改将包含观察到的 y 值的列替换为充满 1 的列（这是与 beta0 相乘所必需的），这就是为什么您没有得到您想要的 2 * 1 矩阵的原因。

其它你可能感兴趣的问题

上一篇如何计算非线性方程的 95% 置信区间？下一篇如果给定一个庞大的数据集，为什么统计模型会过拟合？