连接 r 函数中的多个步骤以计算 beta 值

Question

我是 r 的新手（但我主要用 python 编写代码），我正在尝试写下简单线性回归的代码以供我自己理解，我正处于估计 beta1 的阶段

给定一个预测列 x 和一个响应 y 我想用伪代码做：

sum((x[i] - mean(x)) * (y[i] - mean(y)) / sum(x[i] - mean(x))^2

所以在 r 中：

m <- rbind(c(2,3),c(1,2),c(0,3))) 

自从我读到 for 循环是魔鬼......我想也许我可以做这样的事情：

beta1 <- function(x, y){
    c <- cbind(x,y)
    b1 <- apply(c, 2, function(v) v - mean(v))
    b1 <- b1[,1] * b1[,2]
    b1top <- sum(b1)
    b1bottom <- sum((x - mean(x))^2)
    b1 <- b1top / b1bottom
    return(b1)
}

beta1(m)
[1] 0

现在，抛开实现可能开始错误不谈......有哪些方法可以缩短工作量，就函数内部的代码行而言？

Answer 1

你说得对，for 循环很糟糕。鉴于您尽可能采用矢量化方法，您的方法已经非常快（即您将 x 视为一个向量，并通过从每个元素中减去 x 的平均值来进行标量减法，而不是循环遍历它手动并在循环内减去）。无论如何，这就是你在代码的后半部分所做的事情

您可以用同样的方式缩短代码的前半部分

因此，您可以直接执行此操作，而不是应用该函数来减去平均值（例如x - mean(x)）。这意味着您的分子可以这样计算：

b1Top <- sum((x - mean(x)) * (y - mean(y)))
b1bottom <- sum((x - mean(x))^2)
b1 <- b1top / b1bottom

但是，如果您有超过 1 个预测变量，那么该方法会变得有点繁重。还有另一种使用矢量化方法（使用矩阵）计算回归权重的方法。

可以完全使用原始数据的矩阵运算来计算回归权重。权重由：

X 是您的预测变量矩阵，Y 是您的响应变量

首先，我们需要通过获取预测变量并为截距添加 1 列来创建预测变量/设计矩阵：

xData <- data.frame(1, x)
designMatrix <- data.matrix(xData)

接下来我们计算平方和叉积矩阵（X'X）：

SSCP <- t(designMatrix) %*% designMatrix

然后反转它：

inverseSSCP <- solve(SSCP)

乘以设计矩阵的转置：

inverseMult <- inverseSSCP %*% t(designMatrix)

最后乘以 Y 向量：

betas <- inverseMult %*% y

原始方法不能很好地扩展到多个预测变量，因为您将开始失去 R 的矢量化功能，因此当您到达那个阶段时，您最终将编写更多代码行。矩阵方法允许您一次性计算所有预测变量的回归权重，而不管有多少。