Я хочу вычислить обычные оценки наименьших квадратов ( OLS ) в R без использования «lm» , и это по нескольким причинам.Во-первых, «lm» также вычисляет много вещей, которые мне не нужны (например, подогнанные значения), учитывая, что размер данных является проблемой в моем случае.Во-вторых, я хочу иметь возможность самостоятельно реализовать OLS в R, прежде чем делать это на другом языке (например, в C с GSL).
Как вы, возможно, знаете, модель это:Y = Xb + E;с E ~ N (0, сигма ^ 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 и т. Д. Кроме того, я хочу использовать объясняющие переменные один y один: для выхода 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))
}
Вот мой код, использующий встроенную-in "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)