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