如何在响应矩阵的每一列上计算最小但快速的线性回归？
|
我想在不使用\“ lm \”的情况下计算R中的普通最小二乘（OLS）估计，这有几个原因。首先，考虑到数据大小对于我来说是一个问题，\“ lm \”还计算出许多我不需要的东西（例如拟合值）。其次，我希望能够在使用其他语言（例如，使用GSL的C语言）之前在R中自己实现OLS。
您可能知道，该模型为：Y = Xb + E;与E〜N（0，sigma ^ 2）。如下详述，b是一个具有2个参数的向量，均值（b0）和另一个系数（b1）。最后，对于每次线性回归，我都希望得到b1的估计值（效果大小），标准误差，sigma ^ 2的估计值（残差）和R ^ 2的确定系数。
这是我的数据。我有N个样本（例如，个人，N〜= 100）。对于每个样本，我有Y个输出（响应变量，Y〜= 10 ^ 3）和X个点（解释变量，X〜= 10 ^ 6）。我想分别处理Y输出，即。我要启动Y线性回归：一个用于输出1，一个用于输出2，依此类推。此外，我想使用解释变量1 y 1：对于输出1，我想在点1，然后在点2上对其进行回归。然后...最后在X点。（我希望这很清楚...！）
这是我的R代码，用于检查“ lm”的速度与通过矩阵代数计算OLS估计自己的速度。
首先，我模拟虚拟数据：
nb.samples <- 10 # N
nb.points <- 1000 # X
x <- matrix(data=replicate(nb.samples,sample(x=0:2,size=nb.points, replace=T)),
nrow=nb.points, ncol=nb.samples, byrow=F,
dimnames=list(points=paste(\"p\",seq(1,nb.points),sep=\"\"),
samples=paste(\"s\",seq(1,nb.samples),sep=\"\")))
nb.outputs <- 10 # Y
y <- matrix(data=replicate(nb.outputs,rnorm(nb.samples)),
nrow=nb.samples, ncol=nb.outputs, byrow=T,
dimnames=list(samples=paste(\"s\",seq(1,nb.samples),sep=\"\"),
outputs=paste(\"out\",seq(1,nb.outputs),sep=\"\")))
GetResFromCustomLinReg <- function(Y, xi){ # both Y and xi are N-dim vectors
n <- length(Y)
X <- cbind(rep(1,n), xi) #
p <- 1 # nb of explanatory variables, besides the mean
r <- p + 1 # rank of X: nb of indepdt explanatory variables
inv.XtX <- solve(t(X) %*% X)
beta.hat <- inv.XtX %*% t(X) %*% Y
Y.hat <- X %*% beta.hat
E.hat <- Y - Y.hat
E2.hat <- (t(E.hat) %*% E.hat)
sigma2.hat <- (E2.hat / (n - r))[1,1]
var.covar.beta.hat <- sigma2.hat * inv.XtX
se.beta.hat <- t(t(sqrt(diag(var.covar.beta.hat))))
Y.bar <- mean(Y)
R2 <- 1 - (E2.hat) / (t(Y-Y.bar) %*% (Y-Y.bar))
return(c(beta.hat[2], se.beta.hat[2], sigma2.hat, R2))
}
res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)})
res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)})
> system.time( res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)}) )
user system elapsed
2.526 0.000 2.528
> system.time( res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
user system elapsed
4.561 0.000 4.561
> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)
没有找到相关结果
已邀请:
2 个回复
淘圃跺枯替
函数。在此之前，RcppGSL中也有类似的“ 6”字样，但是您可能需要基于Armadillo的解决方案。在较早的演示文稿中（在带有R的HPC上），我有一些幻灯片显示了速度的提高。 还要注意，帮助页面中的提示是关于“ \ pivoted \”方法比X \'X的反函数更好的提示，后者对退化的模型矩阵可能会很重要。嘘伪
> system.time( res1 <- apply(x, 1, function(xi){lm(y ~ xi)}) ) user system elapsed 2.859 0.005 2.865 > system.time( res2 <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) ) user system elapsed 4.620 0.004 4.626 > system.time( res3 <- apply(x, 1, function(xi){lm.fit(x=cbind(1,xi), y=y)}) ) user system elapsed 0.454 0.004 0.458 > system.time( res4 <- apply(x, 1, function(xi){apply(y, 2, fastLm, x=cbind(1,xi))}) ) user system elapsed 2.279 0.005 2.283 > system.time( res5 <- apply(x, 1, function(xi){apply(y, 2, fastLmPure, cbind(1,xi))}) ) user system elapsed 1.053 0.003 1.056例如，相对于\“ lm \”，我们可能更喜欢使用\“ lm.fit \”，但是如果我们需要R ^ 2，则必须自己计算。同上，我们可能想在\“ lm \”上使用\“ fastLm \”，但是如果要对sigma进行估算，则必须自己计算。用自定义的R函数计算这些东西可能不是很有效（与\“ lm \”所做的比较）。 考虑到所有这些，我暂时将继续使用“ lm”，但是确实Dirk对“ fastLm”的评论是一个很好的建议（这就是我选择他的答案的原因，因为它应该让其他人感兴趣）。