[Rcpp-commits] r3066 - in pkg/RcppArmadillo: . R inst/unitTests man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 15 00:42:24 CEST 2011
Author: dmbates
Date: 2011-06-15 00:42:24 +0200 (Wed, 15 Jun 2011)
New Revision: 3066
Modified:
pkg/RcppArmadillo/DESCRIPTION
pkg/RcppArmadillo/R/fastLm.R
pkg/RcppArmadillo/inst/unitTests/runit.fastLm.R
pkg/RcppArmadillo/man/fastLm.Rd
pkg/RcppArmadillo/src/fastLm.cpp
Log:
Change order of arguments in fastLm.cpp and fastLm.R and documentation. New version.
Modified: pkg/RcppArmadillo/DESCRIPTION
===================================================================
--- pkg/RcppArmadillo/DESCRIPTION 2011-06-14 22:34:50 UTC (rev 3065)
+++ pkg/RcppArmadillo/DESCRIPTION 2011-06-14 22:42:24 UTC (rev 3066)
@@ -1,7 +1,7 @@
Package: RcppArmadillo
Type: Package
Title: Rcpp integration for Armadillo templated linear algebra library
-Version: 0.2.22
+Version: 0.2.23
Date: $Date$
Author: Romain Francois, Dirk Eddelbuettel and Doug Bates
Maintainer: Romain Francois, Dirk Eddelbuettel and Doug Bates <RcppArmadillo-authors at r-enthusiasts.com>
Modified: pkg/RcppArmadillo/R/fastLm.R
===================================================================
--- pkg/RcppArmadillo/R/fastLm.R 2011-06-14 22:34:50 UTC (rev 3065)
+++ pkg/RcppArmadillo/R/fastLm.R 2011-06-14 22:42:24 UTC (rev 3066)
@@ -18,28 +18,28 @@
## You should have received a copy of the GNU General Public License
## along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
-fastLmPure <- function(y, X) {
+fastLmPure <- function(X, y) {
stopifnot(is.matrix(X))
stopifnot(nrow(y)==nrow(X))
- res <- .Call("fastLm", y, X, package="RcppArmadillo")
+ .Call("fastLm", X, y, package = "RcppArmadillo")
}
-fastLm <- function(x, ...) UseMethod("fastLm")
+fastLm <- function(X, ...) UseMethod("fastLm")
-fastLm.default <- function(x, y, ...) {
+fastLm.default <- function(X, y, ...) {
- x <- as.matrix(x)
+ X <- as.matrix(X)
y <- as.numeric(y)
- res <- fastLmPure(y, x)
+ res <- fastLmPure(X, y)
- res$coefficients <- res$coefficient[,1] # force into single-col vector
+ res$coefficients <- as.vector(res$coefficient)
- names(res$coefficients) <- colnames(x)
+ names(res$coefficients) <- colnames(X)
- res$fitted.values <- as.vector(x %*% res$coefficients)
+ res$fitted.values <- as.vector(X %*% res$coefficients)
res$residuals <- y - res$fitted.values
res$call <- match.call()
@@ -65,7 +65,7 @@
# why do I need this here?
rownames(TAB) <- names(object$coefficients)
- colnames(TAB) <- c("Estimate", "StdErr", "t.value", "p.value")
+# colnames(TAB) <- c("Estimate", "StdErr", "t.value", "p.value")
## cf src/stats/R/lm.R and case with no weights and an intercept
f <- object$fitted.values
Modified: pkg/RcppArmadillo/inst/unitTests/runit.fastLm.R
===================================================================
--- pkg/RcppArmadillo/inst/unitTests/runit.fastLm.R 2011-06-14 22:34:50 UTC (rev 3065)
+++ pkg/RcppArmadillo/inst/unitTests/runit.fastLm.R 2011-06-14 22:42:24 UTC (rev 3066)
@@ -25,8 +25,8 @@
test.fastLm <- function() {
data(trees)
flm <- .Call("fastLm",
+ cbind(1, log(trees$Girth)),
log(trees$Volume),
- cbind(1, log(trees$Girth)),
package="RcppArmadillo")
fit <- lm(log(Volume) ~ log(Girth), data=trees)
Modified: pkg/RcppArmadillo/man/fastLm.Rd
===================================================================
--- pkg/RcppArmadillo/man/fastLm.Rd 2011-06-14 22:34:50 UTC (rev 3065)
+++ pkg/RcppArmadillo/man/fastLm.Rd 2011-06-14 22:42:24 UTC (rev 3066)
@@ -10,16 +10,15 @@
function of \code{Armadillo} linear algebra library.
}
\usage{
-fastLmPure(y, X)
+fastLmPure(X, y)
-fastLm(x, \dots)
-\method{fastLm}{default}(x, y, \dots)
+fastLm(X, \dots)
+\method{fastLm}{default}(X, y, \dots)
\method{fastLm}{formula}(formula, data = list(), \dots)
}
\arguments{
\item{y}{a vector containing the explained variable.}
- \item{X}{a data.frame or matrix containing the explanatory variables.}
- \item{x}{a numeric design matrix for the model.}
+ \item{X}{a model matrix.}
\item{formula}{a symbolic description of the model to be fit.}
\item{data}{an optional data frame containing the variables in the model.}
\item{\ldots}{not used}
@@ -28,8 +27,8 @@
Linear models should be estimated using the \code{\link{lm}} function. In
some cases, \code{\link{lm.fit}} may be appropriate.
- The \code{fastLmPure} function provides a reference use case of the \code{GSL}
- library via the wrapper functions in the \pkg{RcppGSL} package.
+ The \code{fastLmPure} function provides a reference use case of the \code{Armadillo}
+ library via the wrapper functions in the \pkg{RcppArmadillo} package.
The \code{fastLm} function provides a more standard implementation of
a linear model fit, offering both a default and a formula interface as
@@ -57,7 +56,7 @@
\code{fastLmPure} returns a list with three components:
\item{coefficients}{a vector of coefficients}
\item{stderr}{a vector of the (estimated) standard errors of the coefficient estimates}
- \item{df}{a scalar denoting the degrees of freedom in the model}
+ \item{df.residual}{a scalar denoting the degrees of freedom in the model}
\code{fastLm} returns a richer object which also includes the
residuals, fitted values and call argument similar to the \code{\link{lm}} or
@@ -71,10 +70,10 @@
}
\examples{
require(datasets)
- data(trees)
+# data(trees) # Unnecessary - the datasets package uses LazyData
## bare-bones direct interface
- flm <- fastLmPure( log(trees$Volume), cbind(1, log(trees$Girth)) )
+ flm <- fastLmPure( cbind(1, log(trees$Girth)), log(trees$Volume) )
print(flm)
## standard R interface for formula or data returning object of class fastLm
Modified: pkg/RcppArmadillo/src/fastLm.cpp
===================================================================
--- pkg/RcppArmadillo/src/fastLm.cpp 2011-06-14 22:34:50 UTC (rev 3065)
+++ pkg/RcppArmadillo/src/fastLm.cpp 2011-06-14 22:42:24 UTC (rev 3066)
@@ -21,23 +21,26 @@
#include <RcppArmadillo.h>
-extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
+extern "C" SEXP fastLm(SEXP Xs, SEXP ys) {
try {
- arma::colvec y = Rcpp::as<arma::colvec>(ys); // direct from SEXP to arma::mat
- arma::mat X = Rcpp::as<arma::mat>(Xs);
- int df = X.n_rows - X.n_cols;
+ Rcpp::NumericVector yr(ys); // creates Rcpp vector from SEXP
+ Rcpp::NumericMatrix Xr(Xs); // creates Rcpp matrix from SEXP
+ int n = Xr.nrow(), k = Xr.ncol();
+ arma::mat X(Xr.begin(), n, k, false); // reuses memory and avoids extra copy
+ arma::colvec y(yr.begin(), yr.size(), false);
+
arma::colvec coef = arma::solve(X, y); // fit model y ~ X
arma::colvec res = y - X*coef; // residuals
- double s2 = std::inner_product(res.begin(), res.end(), res.begin(), 0.0)/df;
+ double s2 = std::inner_product(res.begin(), res.end(), res.begin(), 0.0)/(n - k);
// std.errors of coefficients
arma::colvec std_err = arma::sqrt(s2 * arma::diagvec( arma::pinv(arma::trans(X)*X) ));
return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
Rcpp::Named("stderr") = std_err,
- Rcpp::Named("df") = df
+ Rcpp::Named("df.residual") = n - k
);
} catch( std::exception &ex ) {
More information about the Rcpp-commits
mailing list