rcpp/html/fastLMviaGSL_8r_source.html

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