[Rcpp-commits] r823 - in pkg/RcppArmadillo: . R inst inst/include inst/unitTests man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 2 16:46:15 CET 2010
Author: edd
Date: 2010-03-02 16:46:15 +0100 (Tue, 02 Mar 2010)
New Revision: 823
Added:
pkg/RcppArmadillo/R/fastLm.R
pkg/RcppArmadillo/man/fastLm.Rd
Modified:
pkg/RcppArmadillo/NAMESPACE
pkg/RcppArmadillo/inst/ChangeLog
pkg/RcppArmadillo/inst/include/RcppArmadilloDefines.h.in
pkg/RcppArmadillo/inst/unitTests/runit.RcppArmadillo.R
pkg/RcppArmadillo/src/RcppArmadillo.cpp
Log:
added fastLm C++ function with R wrapper, manual page, and unit test
Modified: pkg/RcppArmadillo/NAMESPACE
===================================================================
--- pkg/RcppArmadillo/NAMESPACE 2010-03-02 15:35:27 UTC (rev 822)
+++ pkg/RcppArmadillo/NAMESPACE 2010-03-02 15:46:15 UTC (rev 823)
@@ -1,2 +1,3 @@
useDynLib( RcppArmadillo)
+export(fastLm) # the only visible function here
Added: pkg/RcppArmadillo/R/fastLm.R
===================================================================
--- pkg/RcppArmadillo/R/fastLm.R (rev 0)
+++ pkg/RcppArmadillo/R/fastLm.R 2010-03-02 15:46:15 UTC (rev 823)
@@ -0,0 +1,27 @@
+
+## fastLm.R: Rcpp/Armadillo implementation of lm()
+##
+## Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
+##
+## This file is part of RcppArmadillo.
+##
+## RcppArmadillo is free software: you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 2 of the License, or
+## (at your option) any later version.
+##
+## RcppArmadillo is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
+
+fastLm <- function(y, X) {
+
+ stopifnot(is.matrix(X))
+ stopifnot(nrow(y)==nrow(X))
+
+ res <- .Call("fastLm", y, X, package="RcppArmadillo")
+}
Modified: pkg/RcppArmadillo/inst/ChangeLog
===================================================================
--- pkg/RcppArmadillo/inst/ChangeLog 2010-03-02 15:35:27 UTC (rev 822)
+++ pkg/RcppArmadillo/inst/ChangeLog 2010-03-02 15:46:15 UTC (rev 823)
@@ -1,3 +1,11 @@
+2010-03-02 Dirk Eddelbuettel <edd at dexter>
+
+ * src/RcppArmadillo.cpp: Added bare-bones 'fastLm' function
+ * R/fastLm.R: Added fastLm R wrapper for C++ fastLm
+ * man/fastLm.Rd: Added manual page
+ * inst/unitTests/runit.RcppArmadillo.R: Added unit test
+ * NAMESPACE: Made fastLm visible
+
2010-03-02 Romain Francois <romain at r-enthusiasts.com>
* inst/include/RcppArmadillo.h: support for Armadillo 0.9.*
Modified: pkg/RcppArmadillo/inst/include/RcppArmadilloDefines.h.in
===================================================================
--- pkg/RcppArmadillo/inst/include/RcppArmadilloDefines.h.in 2010-03-02 15:35:27 UTC (rev 822)
+++ pkg/RcppArmadillo/inst/include/RcppArmadilloDefines.h.in 2010-03-02 15:46:15 UTC (rev 823)
@@ -20,5 +20,5 @@
// You should have received a copy of the GNU General Public License
// along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
-#define ARMA_HAS_CUBE @ARMA_HAS_CUBE@
+#define ARMA_VERSION_GE_070 @ARMA_VERSION_GE_070@
#define ARMA_VERSION_GE_090 @ARMA_VERSION_GE_090@
Modified: pkg/RcppArmadillo/inst/unitTests/runit.RcppArmadillo.R
===================================================================
--- pkg/RcppArmadillo/inst/unitTests/runit.RcppArmadillo.R 2010-03-02 15:35:27 UTC (rev 822)
+++ pkg/RcppArmadillo/inst/unitTests/runit.RcppArmadillo.R 2010-03-02 15:46:15 UTC (rev 823)
@@ -19,21 +19,21 @@
test.wrap.R <- function(){
res <- .Call( "RcppArmadillo_wrap", PACKAGE = "RcppArmadillo" )
-
+
checkEquals( res[[1]][[1]], matrix(as.integer((diag(3))),nr=3), msg = "eye<imat>(3,3)" )
checkEquals( res[[1]][[2]], diag(3), msg = "eye<mat>(3,3)" )
checkEquals( res[[1]][[3]], diag(3), msg = "eye<fmat>(3,3)" )
checkEquals( res[[1]][[4]], matrix(as.integer((diag(3))),nr=3), msg = "eye<umat>(3,3)" )
-
+
checkEquals( res[[2]][[1]], matrix(0, ncol = 5, nrow=1), msg = "zeros<mat>(5,1)" )
checkEquals( res[[2]][[2]], matrix(0, ncol = 5, nrow=1), msg = "zeros<fmat>(5,1)" )
-
+
checkEquals( res[[3]][[1]], matrix(0, ncol = 1, nrow=5), msg = "zeros<mat>(1,5)" )
checkEquals( res[[3]][[2]], matrix(0, ncol = 1, nrow=5), msg = "zeros<mat>(1,5)" )
-
+
checkEquals( res[[4]][[1]], matrix(0:3, ncol = 2, nrow=2), msg = "field<int>" )
checkEquals( res[[4]][[2]], matrix(letters[1:4], ncol = 2, nrow=2), msg = "field<std::string>" )
-}
+}
test.wrap.Glue <- function(){
res <- .Call( "RcppArmadillo_wrap_Glue", PACKAGE = "RcppArmadillo" )
@@ -46,7 +46,7 @@
}
test.as.Mat <- function(){
-
+
integer_mat <- matrix( as.integer(diag(4)), ncol = 4, nrow = 4 )
numeric_mat <- diag(5)
res <- .Call( "RcppArmadillo_as_Mat",
@@ -69,3 +69,17 @@
checkEquals( unlist( res ), rep(55.0, 4 ), msg = "as<Row>" )
}
+test.fastLm <- function() {
+ data(trees)
+ flm <- .Call("fastLm",
+ log(trees$Volume),
+ cbind(rep(1,31), log(treesGirth)),
+ PACKAGE="RcppArmadillo")
+ fit <- lm(log(Volume) ~ log(Girth), data=trees)
+
+ checkEquals( as.numeric(flm$coef), as.numeric(coef(fit)),
+ msg="fastLm.coef")
+ checkEquals( as.numeric(flm$stderr), as.numeric(coef(summary(fit))[,2]),
+ msg="fastLm.stderr")
+}
+
Added: pkg/RcppArmadillo/man/fastLm.Rd
===================================================================
--- pkg/RcppArmadillo/man/fastLm.Rd (rev 0)
+++ pkg/RcppArmadillo/man/fastLm.Rd 2010-03-02 15:46:15 UTC (rev 823)
@@ -0,0 +1,44 @@
+\name{fastLm}
+\alias{fastLm}
+\concept{regression}
+\title{Bare-bones linear model fitting function}
+\description{
+ \code{fastLm} estimates the linear model using the \code{solve}
+ function of \code{Armadillo} linear algebra library.
+}
+\usage{
+fastLm(y, X)
+}
+\arguments{
+ \item{y}{a vector containing the explained variable.}
+ \item{X}{a data.frame or matrix containing the explanatory variables.}
+}
+\details{
+ Linear models should be estimated using the \code{\link{lm}} function. In
+ some cases, \code{\link{lm.fit}} may be appropriate.
+
+ This function provides a reference use case of the \code{Armadillo}
+ library.
+
+ It does not provide a formula interface. It does check for missing
+ values. It does (yet?) return a proper object with suitable accessor
+ functions. That's why we call it bare-bones.
+}
+\value{
+ \code{fastLm} returns a list with two components:
+ \item{coef}{a vector of coefficients}
+ \item{stderr}{a vector of the (estimated_ standard error of the
+ coefficient estimate}
+}
+\seealso{\code{\link{lm}}, \code{\link{lm.fit}}}
+\references{Armadillo project: \url{http://arma.sourceforge.net/}}
+\author{
+ Armadillo is written by Conrad Sanderson. RcppArmadillo is written by
+ Romain Francois and Dirk Eddelbuettel.
+}
+\examples{
+ data(trees)
+ fastLm( log(trees$Volume), cbind(rep(1,31), log(trees$Girth)) )
+}
+\keyword{regression}
+
Modified: pkg/RcppArmadillo/src/RcppArmadillo.cpp
===================================================================
--- pkg/RcppArmadillo/src/RcppArmadillo.cpp 2010-03-02 15:35:27 UTC (rev 822)
+++ pkg/RcppArmadillo/src/RcppArmadillo.cpp 2010-03-02 15:46:15 UTC (rev 823)
@@ -1,4 +1,26 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+//
+// RcppArmadillo.cpp: Rcpp/Armadillo glue
+//
+// Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
+//
+// This file is part of RcppArmadillo.
+//
+// RcppArmadillo is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// RcppArmadillo is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
+
#include <RcppArmadillo.h>
+
using namespace Rcpp ;
using namespace arma ;
@@ -156,3 +178,30 @@
return res ;
}
+
+extern "C" SEXP fastLm(SEXP sy, SEXP sX) {
+
+ Rcpp::NumericVector yr(sy);
+ Rcpp::NumericMatrix Xr(sX);
+ std::vector<int> dims = Xr.attr("dim") ;
+ int n = dims[0], k = dims[1];
+
+#if ARMA_VERSION_GE_090
+ arma::mat X(Xr.begin(), n, k, false); // reuses memory and avoids extra copy
+#else
+ arma::mat X(Xr.begin(), n, k);
+#endif
+
+ arma::colvec y(yr.begin(), yr.size());
+
+ arma::colvec coef = solve(X, y); // fit model y ~ X
+
+ arma::colvec resid = y - X*coef;
+ double sig2 = trans(resid)*resid/(n-k);
+ arma::mat covmat = sig2 * arma::inv(arma::trans(X)*X);
+ arma::colvec stderrest = sqrt(covmat.diag());
+
+ Rcpp::Pairlist res(Rcpp::Named( "coef", coef),
+ Rcpp::Named( "stderr", stderrest));
+ return res;
+}
More information about the Rcpp-commits
mailing list