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