Rcpp Version 0.12.12
fastLMviaGSL.r
Go to the documentation of this file.
1 #!/usr/bin/env r
2 #
3 # A faster lm() replacement based on GNU GSL
4 #
5 # This first appeared in the 'Intro to HPC tutorials'
6 # but has been wrapped in inline::cfunction() here
7 #
8 # Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
9 #
10 # This file is part of Rcpp.
11 #
12 # Rcpp is free software: you can redistribute it and/or modify it
13 # under the terms of the GNU General Public License as published by
14 # the Free Software Foundation, either version 2 of the License, or
15 # (at your option) any later version.
16 #
17 # Rcpp is distributed in the hope that it will be useful, but
18 # WITHOUT ANY WARRANTY; without even the implied warranty of
19 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 # GNU General Public License for more details.
21 #
22 # You should have received a copy of the GNU General Public License
23 # along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
24 
25 source("lmGSL.R")
26 
27 checkLmGSL <- function(y, X) {
28  fun <- lmGSL()
29  res <- fun(y, X)
30  fit <- lm(y ~ X - 1)
31  rc <- all.equal( res[[1]], as.numeric(coef(fit))) &
32  all.equal( res[[2]], as.numeric(coef(summary(fit))[,2]))
33  invisible(rc)
34 }
35 
36 timeLmGSL <- function(y, X, N) {
37  fun <- lmGSL();
38  meantime <- mean(replicate(N, system.time(fun(y, X))["elapsed"]), trim=0.05)
39 }
40 
41 set.seed(42)
42 n <- 5000
43 k <- 9
44 X <- cbind( rep(1,n), matrix(rnorm(n*k), ncol=k) )
45 truecoef <- 1:(k+1)
46 y <- as.numeric(X %*% truecoef + rnorm(n))
47 
48 N <- 100
49 
50 stopifnot(checkLmGSL(y, X))
51 mt <- timeLmGSL(y, X, N)
52 cat("GSL: Running", N, "simulations yields (trimmed) mean time", mt, "\n")