rcpp/html/fastLMviaArmadillo_8r_source.html

00001 #!/usr/bin/r -t
00002 #
00003 # A faster lm() replacement based on Armadillo
00004 #
00005 # This improves on the previous version using GNU GSL
00006 #
00007 # Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
00008 #
00009 # This file is part of Rcpp.
00010 #
00011 # Rcpp is free software: you can redistribute it and/or modify it
00012 # under the terms of the GNU General Public License as published by
00013 # the Free Software Foundation, either version 2 of the License, or
00014 # (at your option) any later version.
00015 #
00016 # Rcpp is distributed in the hope that it will be useful, but
00017 # WITHOUT ANY WARRANTY; without even the implied warranty of
00018 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019 # GNU General Public License for more details.
00020 #
00021 # You should have received a copy of the GNU General Public License
00022 # along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
00023
00024 source("lmArmadillo.R")
00025
00026 checkLmArmadillo <- function(y, X) {
00027     fun <- lmArmadillo()
00028     res <- fun(y, X)
00029     fit <- lm(y ~ X - 1)
00030     rc <- all.equal( res[[1]], as.numeric(coef(fit))) &
00031           all.equal( res[[2]], as.numeric(coef(summary(fit))[,2]))
00032     invisible(rc)
00033 }
00034
00035 timeLmArmadillo <- function(y, X, N) {
00036     fun <- lmArmadillo();
00037     meantime <- mean(replicate(N, system.time(fun(y, X))["elapsed"]), trim=0.05)
00038 }
00039
00040 set.seed(42)
00041 n <- 5000
00042 k <- 9
00043 X <- cbind( rep(1,n), matrix(rnorm(n*k), ncol=k) )
00044 truecoef <- 1:(k+1)
00045 y <- as.numeric(X %*% truecoef + rnorm(n))
00046
00047 N <- 100
00048
00049 stopifnot(checkLmArmadillo(y, X))
00050 mt <- timeLmArmadillo(y, X, N)
00051 cat("Armadillo: Running", N, "simulations yields (trimmed) mean time", mt, "\n")