[Rcpp-commits] r1236 - in pkg/RcppArmadillo: . R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 14 02:58:17 CEST 2010
Author: edd
Date: 2010-05-14 02:58:17 +0200 (Fri, 14 May 2010)
New Revision: 1236
Modified:
pkg/RcppArmadillo/NAMESPACE
pkg/RcppArmadillo/NEWS
pkg/RcppArmadillo/R/fastLm.R
pkg/RcppArmadillo/inst/ChangeLog
pkg/RcppArmadillo/man/fastLm.Rd
pkg/RcppArmadillo/src/fastLm.cpp
Log:
RcppArmadillo gets and expandedn fastLm formula interface and summary/print/preditc methods
Modified: pkg/RcppArmadillo/NAMESPACE
===================================================================
--- pkg/RcppArmadillo/NAMESPACE 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/NAMESPACE 2010-05-14 00:58:17 UTC (rev 1236)
@@ -1,5 +1,11 @@
useDynLib(RcppArmadillo)
-export(fastLm)
+export(fastLmCall,
+ fastLm,
+ fastLm.default,
+ fastLm.formula,
+ print.fastLm,
+ summary.fastLm,
+ print.summary.fastLm)
export(RcppArmadillo.package.skeleton)
Modified: pkg/RcppArmadillo/NEWS
===================================================================
--- pkg/RcppArmadillo/NEWS 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/NEWS 2010-05-14 00:58:17 UTC (rev 1236)
@@ -1,19 +1,23 @@
0.2.0 (under development)
+ o fastLm is now generic and has a formula interface as well as methods
+ for print, summary, predict to behave like a standard model fitting
+ function
+
o Armadillo sources (using release 0.9.6) are now included in the package
- using a standardized build suitable for our purposes (not assuming
- Boost or Atlas) -- see ?RcppArmadillo for details
+ using a standardized build suitable for our purposes (not assuming
+ Boost or Atlas) -- see ?RcppArmadillo for details
o New R function RcppArmadillo.package.skeleton, similar to
- Rcpp::Rcpp.package.skeleton, but targetting use of RcppArmadillo
+ Rcpp::Rcpp.package.skeleton, but targetting use of RcppArmadillo
0.1.0 2010-03-11
o the fastLm() implementation of a bare-bones lm() fit
- (using Armadillo's solve() function) provides an example
- of how efficient code can be written compactly using
- the combination of Rcpp, RcppAramadillo and Armadillo
+ (using Armadillo's solve() function) provides an example
+ of how efficient code can be written compactly using
+ the combination of Rcpp, RcppAramadillo and Armadillo
o support for Rcpp implicit wrap of these types :
Mat<T>, Col<T>, Row<T>, Cube<T> where T is one of :
@@ -23,3 +27,4 @@
Mat<T>, Col<T>, Row<T> where R is one of :
int, unsigned int, double, float
+
Modified: pkg/RcppArmadillo/R/fastLm.R
===================================================================
--- pkg/RcppArmadillo/R/fastLm.R 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/R/fastLm.R 2010-05-14 00:58:17 UTC (rev 1236)
@@ -18,7 +18,7 @@
## 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) {
+fastLmCall <- function(y, X) {
stopifnot(is.matrix(X))
stopifnot(nrow(y)==nrow(X))
@@ -26,18 +26,82 @@
res <- .Call("fastLm", y, X, package="RcppArmadillo")
}
-## What would be nice here:
-##
-## fastLm <- function(x, ...) UseMethod("fastLm")
-##
-## fastLm.formula <- ...
-##
-## fastLm.default <- ...
-##
-## print.fastLm <- function(x, ...)
-##
-## summary.fastLm <- function(object, ...
-##
-## print.summary.fastLm <- ...
-##
-## but on the other hand lm.fit() does not do any of this either
+fastLm <- function(x, ...) UseMethod("fastLm")
+
+fastLm.default <- function(x, y, ...) {
+
+ x <- as.matrix(x)
+ y <- as.numeric(y)
+
+ res <- fastLmCall(y, x)
+
+ names(res$coefficients) <- colnames(res$vcov) <- rownames(res$vcov) <- colnames(x)
+
+ res$fitted.values <- as.vector(x %*% res$coefficients)
+ res$residuals <- y - res$fitted.values
+ res$call <- match.call()
+
+ class(res) <- "fastLm"
+ res
+}
+
+print.fastLm <- function(x, ...) {
+ cat("\nCall:\n")
+ print(x$call)
+ cat("\nCoefficients:\n")
+ print(x$coefficients)
+}
+
+summary.fastLm <- function(object, ...) {
+ se <- object$stderr
+ tval <- coef(object)/se
+
+ TAB <- cbind(Estimate = coef(object),
+ StdErr = se,
+ t.value = tval,
+ p.value = 2*pt(-abs(tval), df=object$df))
+
+ # why do I need this here?
+ rownames(TAB) <- names(object$coefficients)
+ colnames(TAB) <- c("Estimate", "StdErr", "t.value", "p.value")
+
+ res <- list(call=object$call,
+ coefficients=TAB)
+
+ class(res) <- "summary.fastLm"
+ res
+}
+
+print.summary.fastLm <- function(x, ...) {
+ cat("\nCall:\n")
+ print(x$call)
+ cat("\n")
+
+ printCoefmat(x$coefficients, P.value=TRUE, has.Pvalue=TRUE)
+}
+
+fastLm.formula <- function(formula, data=list(), ...) {
+ mf <- model.frame(formula=formula, data=data)
+ x <- model.matrix(attr(mf, "terms"), data=mf)
+ y <- model.response(mf)
+
+ res <- fastLm.default(x, y, ...)
+ names(res$coefficients) <- colnames(x) <- colnames(res$vcov) <- rownames(res$vcov) <- colnames(x)
+ res$call <- match.call()
+ res$formula <- formula
+ res
+}
+
+predict.fastLm <- function(object, newdata=NULL, ...) {
+ if (is.null(newdata)) {
+ y <- fitted(object)
+ } else {
+ if (!is.null(object$formula)) {
+ x <- model.matrix(object$formula, newdata)
+ } else {
+ x <- newdata
+ }
+ y <- as.vector(x %*% coef(object))
+ }
+ y
+}
Modified: pkg/RcppArmadillo/inst/ChangeLog
===================================================================
--- pkg/RcppArmadillo/inst/ChangeLog 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/inst/ChangeLog 2010-05-14 00:58:17 UTC (rev 1236)
@@ -1,3 +1,14 @@
+2010-05-13 Dirk Eddelbuettel <edd at dexter>
+
+ * R/fastLm.R: fastLm() is now generic and has a formula interface as
+ well as methods for print, summary, predict to behave like a standard
+ model fitting function
+
+ * R/fastLm.Rd: Added documentation for new methods
+
+ * src/fastLm.cpp: Small additions to returned parameters to support
+ new functionality
+
2010-05-04 Romain Francois <romain at r-enthusiasts.com>
* R/RcppArmadillo.package.skeleton.R: new package skeleton function
@@ -3,5 +14,5 @@
similar to Rcpp::Rcpp.package.skeleton but targetting
use of RcppArmadillo
-
+
* README: removed and promoted to additional documentation in ?RcppArmadillo
Modified: pkg/RcppArmadillo/man/fastLm.Rd
===================================================================
--- pkg/RcppArmadillo/man/fastLm.Rd 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/man/fastLm.Rd 2010-05-14 00:58:17 UTC (rev 1236)
@@ -1,5 +1,12 @@
\name{fastLm}
\alias{fastLm}
+\alias{fastLmCall}
+\alias{fastLm.default}
+\alias{fastLm.formula}
+\alias{print.fastLm}
+\alias{summary.fastLm}
+\alias{print.summary.fastLm}
+\alias{predict.fastLm}
\concept{regression}
\title{Bare-bones linear model fitting function}
\description{
@@ -7,11 +14,25 @@
function of \code{Armadillo} linear algebra library.
}
\usage{
-fastLm(y, X)
+fastLmCall(y, X)
+
+fastLm(x, ...)
+\method{fastLm}{default}(x, y, ...)
+\method{fastLm}{formula}(formula, data = list(), ...)
+\method{print}{fastLm}(x, ...)
+\method{summary}{fastLm}(object, ...)
+\method{predict}{fastLm}(object, newdata=NULL, ...)
}
\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{formula}{a symbolic description of the model to be fit.}
+ \item{data}{an optional data frame containing the variables in the model.}
+ \item{newdata}{an optional data frame containing the variables for
+ which a model prediction is desired.}
+ \item{object}{an object of class \code{"fastLm"}, i.e. a fitted model}
+ \item{\ldots}{not used}
}
\details{
Linear models should be estimated using the \code{\link{lm}} function. In
@@ -55,8 +76,19 @@
Romain Francois, Dirk Eddelbuettel and Douglas Bates.
}
\examples{
+ require(datasets)
data(trees)
- flm <- fastLm( log(trees$Volume), cbind(rep(1,31), log(trees$Girth)) )
+
+ ## bare-bones direct interface
+ flm <- fastLmCall( log(trees$Volume), cbind(1, log(trees$Girth)) )
+ print(flm)
+
+ ## standard R interface for formula or data returning object of class fastLm
+ flmmod <- fastLm( log(Volume) ~ log(Girth), data=trees)
+ summary(flmmod)
}
\keyword{regression}
+
+
+
Modified: pkg/RcppArmadillo/src/fastLm.cpp
===================================================================
--- pkg/RcppArmadillo/src/fastLm.cpp 2010-05-13 19:58:26 UTC (rev 1235)
+++ pkg/RcppArmadillo/src/fastLm.cpp 2010-05-14 00:58:17 UTC (rev 1236)
@@ -23,24 +23,36 @@
extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
- 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();
+ try {
+ 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::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 coef = arma::solve(X, y); // fit model y ~ X
arma::colvec resid = y - X*coef; // residuals
double sig2 = arma::as_scalar( arma::trans(resid)*resid/(n-k) );
- // std.error of estimate
- arma::colvec stderrest = arma::sqrt( sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
+
+ arma::mat covmat = sig2 * arma::inv(arma::trans(X)*X); // covmat
+ arma::colvec stderrest = arma::sqrt(arma::diagvec(covmat)); // std.error of estimate
- return Rcpp::List::create(
- Rcpp::Named("coefficients") = coef,
- Rcpp::Named("stderr") = stderrest
- ) ;
-
+ return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
+ Rcpp::Named("stderr") = stderrest,
+ Rcpp::Named("vcov") = covmat,
+ Rcpp::Named("df") = n - k
+ );
+
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
}
+
+
+
+
More information about the Rcpp-commits
mailing list