我想计算普通最小二乘（OLS) R 中的估计不使用“lm”，这有几个原因。首先，考虑到数据大小在我的情况下是一个问题，“lm”还计算了很多我不需要的东西（例如拟合值）。其次，我希望能够在使用另一种语言（例如使用 GSL 的 C 语言）之前先在 R 中实现 OLS。
如您所知，model为：Y=Xb+E；与 E ~ N(0, sigma^2)。如下所述，b 是一个具有 2 个参数的向量，即均值 (b0) 和另一个系数 (b1)。最后，对于我要做的每个线性回归，我想要 b1 （效应大小）的估计值、其标准误差、sigma^2 （残差方差）和 R^2 （确定系数）的估计值。
这是我的data。我有 N 个样本（例如个体，N~=100）。对于每个样本，我有 Y 输出（响应变量，Y~=10^3）和 X 点（解释变量，X~=10^6）。我想单独处理 Y 输出，即。我想启动 Y 线性回归：一个用于输出 1，一个用于输出 2，等等。此外，我想使用解释变量 1 y 1：对于输出 1，我想在点 1 上对其进行回归，然后在点 2 上进行回归，然后......最后关于X点。（我希望它很清楚......！）
这是我的R code检查“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))
}
这是我使用内置“lm”的代码：
res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)})
这是我的自定义 OLS 代码：
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
（当然，当增加 N、X 和 Y 时，情况会变得更糟。）
当然，“lm”具有“自动”单独拟合响应矩阵（y〜xi）的每一列的良好特性，而我必须使用“应用”Y次（对于每个yi〜xi）。但这是我的代码速度较慢的唯一原因吗？你们中有人知道如何改进这一点吗？
（很抱歉问了这么长的问题，但我真的试图提供一个最小但全面的示例。）
> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)
