[Rcpp-commits] r2635 - in pkg/RcppModels: . R inst inst/unitTests man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 1 03:08:34 CET 2010
Author: dmbates
Date: 2010-12-01 03:08:30 +0100 (Wed, 01 Dec 2010)
New Revision: 2635
Added:
pkg/RcppModels/ChangeLog
pkg/RcppModels/R/fastGlm.R
pkg/RcppModels/R/zzz.R
pkg/RcppModels/inst/
pkg/RcppModels/inst/unitTests/
pkg/RcppModels/inst/unitTests/runTests.R
pkg/RcppModels/inst/unitTests/runit.link.R
pkg/RcppModels/tests/
pkg/RcppModels/tests/doRUnit.R
pkg/RcppModels/tests/pois.R
Modified:
pkg/RcppModels/DESCRIPTION
pkg/RcppModels/NAMESPACE
pkg/RcppModels/man/RccpModels-package.Rd
pkg/RcppModels/src/Module.cpp
pkg/RcppModels/src/RcppModels.cpp
pkg/RcppModels/src/RcppModels.h
pkg/RcppModels/src/glmFamily.cpp
pkg/RcppModels/src/glmFamily.h
Log:
see ChangeLog
Added: pkg/RcppModels/ChangeLog
===================================================================
--- pkg/RcppModels/ChangeLog (rev 0)
+++ pkg/RcppModels/ChangeLog 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,27 @@
+2010-11-30 Douglas Bates <bates at stat.wisc.edu>
+
+ * src/Module.cpp: definition of the Rcpp module
+
+ * src/glmFamily.cpp, src/glmFamily.h: A class representing a glm
+ family.
+
+ * src/RcppModels.cpp, src/RcppModels.h: Classes for linear and
+ generalized linear response modules and for dense predictor
+ modules.
+
+ * tests/pois.R: An example from Dobson. The Gauss-Newton version
+ of IRLS is currently being done by hand. This will change.
+
+ * man/RccpModels-package.Rd: initial package documentation
+
+ * inst/unitTests: installed unit tests of some of the glmFamily
+ methods.
+
+ * R/fastGlm.R (fastGlm): draft glm function using C++Class
+ instances. Have not yet added the Gauss-Newton steps but do have
+ one fixed-point step.
+
+ * R/zzz.R (.onLoad): populate the namespace with module contents
+
+ * DESCRIPTION, NAMESPACE: initial release
+
Modified: pkg/RcppModels/DESCRIPTION
===================================================================
--- pkg/RcppModels/DESCRIPTION 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/DESCRIPTION 2010-12-01 02:08:30 UTC (rev 2635)
@@ -1,13 +1,14 @@
Package: RcppModels
Type: Package
Title: Classes for linear and generalized linear and nonlinear models
-Version: 1.0
-Date: 2010-11-20
+Version: 0.1.0
+Date: 2010-11-30
Author: Douglas Bates, John Chambers, Dirk Eddelbuettel, Romain Francois
Maintainer: Rcpp-Core <rcpp-core at r-forge.wu-wien.ac.at>
Description: Provides classes and methods for linear, generalized
linear and nonlinear models that use linear predictor expressions.
License: GPL (>= 2)
LazyLoad: yes
-Depends: RcppArmadillo, Rcpp
+Depends: R(>= 2.12.0), Rcpp(>= 0.8.9), RcppArmadillo, methods,
LinkingTo: RcppArmadillo, Rcpp
+Suggests: RUnit
Modified: pkg/RcppModels/NAMESPACE
===================================================================
--- pkg/RcppModels/NAMESPACE 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/NAMESPACE 2010-12-01 02:08:30 UTC (rev 2635)
@@ -1,3 +1,13 @@
useDynLib(RcppModels)
+import(Rcpp)
+export(
+ densePred
+ , fastGlm
+ , glmFamily
+ , glmResp
+ , lmResp
+ )
+# exportPattern("^[^\\.]")
+
Added: pkg/RcppModels/R/fastGlm.R
===================================================================
--- pkg/RcppModels/R/fastGlm.R (rev 0)
+++ pkg/RcppModels/R/fastGlm.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,67 @@
+mkRespMod <- function(fr, family)
+{
+ N <- n <- nrow(fr)
+ # components of the model frame
+ y <- model.response(fr)
+ if(length(dim(y)) == 1) { # avoid problems with 1D arrays, but keep names
+ nm <- rownames(y)
+ dim(y) <- NULL
+ if(!is.null(nm)) names(y) <- nm
+ }
+ weights <- unname(model.weights(fr))
+ if (is.null(weights)) weights <- rep.int(1, n)
+ offset <- unname(model.offset(fr))
+ if (is.null(offset)) offset <- rep(0, n)
+ rho <- new.env()
+ rho$etastart <- model.extract(fr, "etastart")
+ rho$mustart <- model.extract(fr, "mustart")
+ rho$nobs <- n
+ rho$y <- y
+ if (!"family" %in% class(family))
+ stop("family should be a list with class \"family\"")
+ eval(family$initialize, rho)
+ family$initialize <- NULL # remove clutter from str output
+ rr <- new(glmResp, family, rho$y, weights, offset, rho$n)
+ rr$updateMu(family$linkfun(rho$mustart))
+ rr$updateWts()
+ rr
+}
+
+fastGlm <-
+ function(formula, family, data, weights, subset,
+ na.action, start = NULL, etastart, mustart, offset,
+ drop.unused.levels = FALSE, doFit = TRUE,
+ control = list(...))
+{
+ call <- match.call()
+ if (missing(family)) {
+ family <- NULL
+ } else {
+ if(is.character(family))
+ family <- get(family, mode = "function", envir = parent.frame())
+ if(is.function(family)) family <- family()
+ if(is.null(family$family)) {
+ print(family)
+ stop("'family' not recognized")
+ }
+ }
+ ## extract x, y, etc from the model formula and frame
+ if(missing(data)) data <- environment(formula)
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data", "subset", "weights", "na.action",
+ "etastart", "mustart", "offset"), names(mf), 0L)
+ mf <- mf[c(1L, m)]
+ mf$drop.unused.levels <- TRUE
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+
+ rr <- mkRespMod(mf, family)
+ pp <- new(densePred, model.matrix(formula, mf))
+ ## one iteration of the fixed-point algorithm to establish
+ ## a baseline coefficient vector
+ wts <- rr$sqrtWrkWt
+ rr$updateMu(pp$gamma(wts, wts * rr$wrkResp))
+ rr$updateWts()
+ pp$installCoef0()
+ list(rmod = rr, pmod = pp)
+}
Added: pkg/RcppModels/R/zzz.R
===================================================================
--- pkg/RcppModels/R/zzz.R (rev 0)
+++ pkg/RcppModels/R/zzz.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,6 @@
+NAMESPACE <- environment()
+.onLoad <- function(libname, pkgname) {
+ ## load the module and store it in our namespace
+ mod <- Module("RcppModels")
+ populate(mod, NAMESPACE)
+}
Added: pkg/RcppModels/inst/unitTests/runTests.R
===================================================================
--- pkg/RcppModels/inst/unitTests/runTests.R (rev 0)
+++ pkg/RcppModels/inst/unitTests/runTests.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,88 @@
+pkg <- "RcppModels"
+
+if(require("RUnit", quietly = TRUE)) {
+
+ is_local <- function(){
+ if( exists( "argv", globalenv() ) && "--local" %in% argv ) return(TRUE)
+ if( "--local" %in% commandArgs(TRUE) ) return(TRUE)
+ FALSE
+ }
+ if( is_local() ) path <- getwd()
+
+ library(package=pkg, character.only = TRUE)
+ if(!(exists("path") && file.exists(path)))
+ path <- system.file("unitTests", package = pkg)
+
+ ## --- Testing ---
+
+ ## Define tests
+ testSuite <- defineTestSuite(name=paste(pkg, "unit testing"), dirs = path)
+
+ if(interactive()) {
+ cat("Now have RUnit Test Suite 'testSuite' for package '", pkg,
+ "' :\n", sep='')
+ str(testSuite)
+ cat('', "Consider doing",
+ "\t tests <- runTestSuite(testSuite)", "\nand later",
+ "\t printTextProtocol(tests)", '', sep="\n")
+ } else { ## run from shell / Rscript / R CMD Batch / ...
+ ## Run
+ tests <- runTestSuite(testSuite)
+
+ output <- NULL
+
+ process_args <- function(argv){
+ if( !is.null(argv) && length(argv) > 0 ){
+ rx <- "^--output=(.*)$"
+ g <- grep( rx, argv, value = TRUE )
+ if( length(g) ){
+ sub( rx, "\\1", g[1L] )
+ }
+ }
+ }
+
+ # give a chance to the user to customize where he/she wants
+ # the unit tests results to be stored with the --output= command
+ # line argument
+ if( exists( "argv", globalenv() ) ){
+ # littler
+ output <- process_args(argv)
+ } else {
+ # Rscript
+ output <- process_args(commandArgs(TRUE))
+ }
+
+ # if it did not work, try to use /tmp
+ if( is.null(output) ){
+ if( file.exists( "/tmp" ) ){
+ output <- "/tmp"
+ } else{
+ output <- getwd()
+ }
+ }
+
+ ## Print results
+ output.txt <- file.path( output, sprintf("%s-unitTests.txt", pkg))
+ output.html <- file.path( output, sprintf("%s-unitTests.html", pkg))
+
+ printTextProtocol(tests, fileName=output.txt)
+ message( sprintf( "saving txt unit test report to '%s'", output.txt ) )
+
+ ## Print HTML version to a file
+ ## printHTMLProtocol has problems on Mac OS X
+ if (Sys.info()["sysname"] != "Darwin"){
+ message( sprintf( "saving html unit test report to '%s'", output.html ) )
+ printHTMLProtocol(tests, fileName=output.html)
+ }
+
+ ## stop() if there are any failures i.e. FALSE to unit test.
+ ## This will cause R CMD check to return error and stop
+ if(getErrors(tests)$nFail > 0) {
+ stop("one of the unit tests failed")
+ }
+ }
+} else {
+ cat("R package 'RUnit' cannot be loaded -- no unit tests run\n",
+ "for package", pkg,"\n")
+}
+
Added: pkg/RcppModels/inst/unitTests/runit.link.R
===================================================================
--- pkg/RcppModels/inst/unitTests/runit.link.R (rev 0)
+++ pkg/RcppModels/inst/unitTests/runit.link.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,71 @@
+eps <- .Machine$double.eps
+oneMeps <- 1 - eps
+set.seed(1)
+etas <-
+ lapply(list(-8:8, # equal spacing to asymptotic area
+ runif(20, -8, 8), # random sample from wide uniform dist
+ rnorm(20, 0, 8), # random sample from wide normal dist
+ c(-10^30, rnorm(10, 0, 4), 10^30)), as.numeric)
+etapos <-
+ lapply(list(1:20,
+ rexp(20),
+ rgamma(20, 3),
+ pmax(eps, rnorm(20, 2, 1))), as.numeric)
+
+mubinom <-
+ lapply(list(runif(100, 0, 1),
+ rbeta(100, 1, 3),
+ pmin(pmax(eps, rbeta(100, 0.1, 3)), oneMeps),
+ pmin(pmax(eps, rbeta(100, 3, 0.1)), oneMeps)), as.numeric)
+
+tst.lnki <- function(fam, lst) {
+ ff <- glmFamily$new(fam)
+ unlist(lapply(lst, function(x)
+ checkEquals(fam$linkinv(x), ff$linkInv(x))))
+}
+tst.link <- function(fam, lst) {
+ ff <- glmFamily$new(fam)
+ unlist(lapply(lst, function(x)
+ checkEquals(fam$linkfun(x), ff$linkFun(x))))
+}
+tst.muEta <- function(fam, lst) {
+ ff <- glmFamily$new(fam)
+ unlist(lapply(lst, function(x)
+ checkEquals(fam$mu.eta(x), ff$muEta(x))))
+}
+
+test.uncons.lnki <- function() { # check linkinv for unconstrained eta
+ tst.lnki(binomial(), etas) # binomial with default, logit link
+ tst.muEta(binomial(), etas)
+ tst.lnki(binomial("probit"), etas) # binomial with probit link
+ tst.muEta(binomial("probit"), etas)
+ tst.lnki(binomial("cloglog"), etas) # binomial with cloglog link
+ tst.muEta(binomial("cloglog"), etas)
+ tst.lnki(poisson(), etas) # Poisson with default, log link
+ tst.muEta(poisson(), etas)
+ tst.lnki(gaussian(), etas) # Gaussian with default, identity link
+ tst.muEta(gaussian(), etas)
+}
+
+test.pos.lnki <- function() { # check linkinv for positive eta only
+ set.seed(1)
+ tst.lnki(Gamma(), etapos) # gamma family
+ tst.muEta(Gamma(), etapos)
+ tst.lnki(inverse.gaussian(), etapos) # inverse Gaussian
+ tst.muEta(inverse.gaussian(), etapos)
+}
+
+test.binom.link <- function() { # check link for binomial mu
+ tst.link(binomial(), mubinom)
+ tst.link(binomial("probit"), mubinom)
+}
+
+test.pos.link <- function() { # check link for positive mu
+ tst.link(poisson(), etapos)
+ tst.link(Gamma(), etapos)
+ tst.link(inverse.gaussian(), etapos)
+}
+
+test.uncons.link <- function() { # check link for unconstrained mu
+ tst.link(gaussian(), etas)
+}
Modified: pkg/RcppModels/man/RccpModels-package.Rd
===================================================================
--- pkg/RcppModels/man/RccpModels-package.Rd 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/man/RccpModels-package.Rd 2010-12-01 02:08:30 UTC (rev 2635)
@@ -11,12 +11,12 @@
}
\details{
\tabular{ll}{
-Package: \tab RcppModels\cr
-Type: \tab Package\cr
-Version: \tab 1.0\cr
-Date: \tab 2010-11-20\cr
-License: \tab GPL (>= 2)
-LazyLoad: \tab yes\cr
+ Package: \tab RcppModels\cr
+ Type: \tab Package\cr
+ Version: \tab 0.1.0\cr
+ Date: \tab 2010-11-30\cr
+ License: \tab GPL (>= 2)\cr
+ LazyLoad: \tab yes\cr
}
}
\author{
Modified: pkg/RcppModels/src/Module.cpp
===================================================================
--- pkg/RcppModels/src/Module.cpp 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/src/Module.cpp 2010-12-01 02:08:30 UTC (rev 2635)
@@ -5,24 +5,23 @@
RCPP_MODULE(RcppModels) {
-class_<Irwls>( "Irwls" )
+// class_<Irwls>( "Irwls" )
- .constructor()
- .constructor<NumericMatrix,NumericVector>()
+// .constructor<NumericMatrix,NumericVector>()
- .field_readonly("x", &Irwls::x)
- .field_readonly("y", &Irwls::y)
+// .field_readonly("x", &Irwls::x)
+// .field_readonly("y", &Irwls::y)
- .method("fit", &Irwls::fit)
- ;
+// .method("fit", &Irwls::fit)
+// ;
using namespace glm;
class_<glmFamily>("glmFamily")
- .constructor()
.constructor<List>()
+ .property("family", &glmFamily::fam)
.method("linkFun", &glmFamily::linkFun)
.method("linkInv", &glmFamily::linkInv)
.method("muEta", &glmFamily::muEta)
@@ -32,7 +31,6 @@
class_<lmResp>("lmResp")
- .constructor()
.constructor<NumericVector>()
.constructor<NumericVector,NumericVector>()
.constructor<NumericVector,NumericVector,NumericVector>()
@@ -52,52 +50,75 @@
class_<glmResp>("glmResp")
- .constructor()
+// .constructor()
.constructor<List,NumericVector>()
.constructor<List,NumericVector,NumericVector>()
.constructor<List,NumericVector,NumericVector,NumericVector>()
.constructor<List,NumericVector,NumericVector,NumericVector,NumericVector>()
- .property("eta", &glmResp::eta)
- .property("mu", &glmResp::mu)
- .property("offset", &glmResp::offset)
- .property("sqrtXwt", &glmResp::sqrtXwt)
- .property("sqrtrwt", &glmResp::sqrtrwt)
- .property("weights", &glmResp::weights)
- .property("wtres", &glmResp::wtres)
- .property("y", &glmResp::y)
- .property("wrss", &glmResp::wrss)
- .property("wrkResids", &glmResp::wrkResids)
+ .property("devResid", &glmResp::devResid,
+ "deviance residuals")
+ .property("eta", &glmResp::eta,
+ "current value of the linear predictor")
+ .property("family", &glmResp::family,
+ "name of the glm family")
+ .property("link", &glmResp::link,
+ "name of the glm link")
+ .property("mu", &glmResp::mu,
+ "current value of the mean vector")
+ .property("muEta", &glmResp::muEta,
+ "current value of d mu/d eta")
+ .property("offset", &glmResp::offset,
+ "offset vector - const")
+ .property("residDeviance", &glmResp::residDeviance,
+ "sum of the deviance residuals")
+ .property("sqrtXwt", &glmResp::sqrtXwt,
+ "square root of the weights applied to the model matrix")
+ .property("sqrtrwt", &glmResp::sqrtrwt,
+ "square root of the weights applied to the residuals")
+ .property("sqrtWrkWt", &glmResp::sqrtrwt,
+ "square root of the weights applied to the working response")
+ .property("weights", &glmResp::weights,
+ "prior weights - const")
+ .property("variance", &glmResp::variance,
+ "current (unscaled) variances")
+ .property("wtres", &glmResp::wtres,
+ "current value of the weighted residuals")
+ .property("wrss", &glmResp::wrss,
+ "weighted residual sum of squares")
+ .property("wrkResids", &glmResp::wrkResids,
+ "working residuals - on the eta scale")
+ .property("wrkResp", &glmResp::wrkResp,
+ "working response - on the eta scale")
+ .property("y", &glmResp::y,
+ "numeric response vector - const")
- .method("devResid", &glmResp::devResid)
- .method("residDeviance", &glmResp::residDeviance)
- .method("updateMu", &glmResp::updateMu)
- .method("updateWts", &glmResp::updateWts)
- .method("wrkResids", &glmResp::wrkResids)
+ .method("updateMu", &glmResp::updateMu,
+ "update mu and derived quantities from a new value of eta")
+ .method("updateWts", &glmResp::updateWts,
+ "update the residual and X weights.")
;
-class_<predMod>("predMod")
-
- .constructor()
- .constructor<int>()
- .constructor<NumericVector>()
-
-// For some reason this doesn't work.
-// .property("coef0", &predMod::coef0, &predMod::setCoef0)
- .property("coef0", &predMod::coef0)
- ;
-
class_<densePred>("densePred")
- .constructor()
.constructor<NumericMatrix>()
.constructor<NumericMatrix, NumericVector>()
- .property("coef0", &densePred::coef0)
- .property("delta", &densePred::delta)
- .property("X", &densePred::X)
+ .property("coef", &densePred::coef,
+ "last coefficient vector evaluated")
+ .property("coef0", &densePred::coef0, &densePred::setCoef0,
+ "base coefficient vector")
+ .property("delta", &densePred::delta,
+ "increment in coefficient vector")
+ .property("X", &densePred::X,
+ "dense model matrix")
- .method("gamma", &densePred::gamma)
+ .method("gamma", &densePred::gamma,
+ "evaluate delta from Xwts and wtres; return eta for full step factor")
+ .method("dgamma", &densePred::dgamma,
+ "return eta for specified step factor")
+ .method("installCoef0", &densePred::installCoef0,
+ "install the current coef as coef0 for the next step")
;
}
Modified: pkg/RcppModels/src/RcppModels.cpp
===================================================================
--- pkg/RcppModels/src/RcppModels.cpp 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/src/RcppModels.cpp 2010-12-01 02:08:30 UTC (rev 2635)
@@ -2,23 +2,23 @@
using namespace Rcpp;
-Irwls::Irwls(Rcpp::NumericMatrix xR, Rcpp::NumericVector yR) :
- x(xR), y(yR) {
- int n = x.nrow(), p = x.ncol();
- if (y.size() != n)
- throw std::runtime_error("incompatible dimensions");
- X = arma::mat(x.begin(), n, p, false, true);
- Y = arma::vec(y.begin(), n, false, true);
-}
+// Irwls::Irwls(Rcpp::NumericMatrix xR, Rcpp::NumericVector yR) :
+// x(xR), y(yR) {
+// int n = x.nrow(), p = x.ncol();
+// if (y.size() != n)
+// throw std::runtime_error("incompatible dimensions");
+// X = arma::mat(x.begin(), n, p, false, true);
+// Y = arma::vec(y.begin(), n, false, true);
+// }
-arma::vec Irwls::fit(Rcpp::NumericVector wR) {
- int n = x.nrow();
- if (wR.size() != n)
- throw std::runtime_error("length(weights) != nrow(X)");
- wrt = sqrt(wR);
- arma::vec W(wrt.begin(), n, false, true);
- return arma::solve(diagmat(W) * X, W % Y);
-}
+// arma::vec Irwls::fit(Rcpp::NumericVector wR) {
+// int n = x.nrow();
+// if (wR.size() != n)
+// throw std::runtime_error("length(weights) != nrow(X)");
+// wrt = sqrt(wR);
+// arma::vec W(wrt.begin(), n, false, true);
+// return arma::solve(diagmat(W) * X, W % Y);
+// }
namespace glm {
lmResp::lmResp(Rcpp::NumericVector y)
@@ -116,6 +116,15 @@
return d_fam.devResid(d_mu, d_weights, d_y);
}
+ Rcpp::NumericVector glmResp::muEta() const {
+ return d_fam.muEta(d_eta);
+ }
+
+ Rcpp::NumericVector glmResp::variance() const {
+ return d_fam.variance(d_mu);
+ }
+
+
double glmResp::residDeviance() const {
return sum(devResid());
}
@@ -129,10 +138,18 @@
}
Rcpp::NumericVector glmResp::wrkResids() const {
- // This needs to be modified
- return d_wtres;
+ return (d_y - d_mu)/muEta();
}
+ Rcpp::NumericVector glmResp::wrkResp() const {
+ return d_eta + wrkResids();
+ }
+
+ Rcpp::NumericVector glmResp::sqrtWrkWt() const {
+ Rcpp::NumericVector me = muEta();
+ return sqrt(d_weights * me * me / variance());
+ }
+
double glmResp::updateMu(const Rcpp::NumericVector& gamma) {
d_eta = d_offset + gamma;
NumericVector mmu = d_fam.linkInv(d_eta);
@@ -141,17 +158,30 @@
}
predMod::predMod(int p)
- : d_coef0(p), d_delta(p),
- a_coef0(d_coef0.begin(), d_coef0.size(), false, true),
- a_delta(d_delta.begin(), d_delta.size(), false, true) {
+ : d_coef(p), d_coef0(p), d_delta(p),
+ a_coef (d_coef.begin(), p, false, true),
+ a_coef0(d_coef0.begin(), p, false, true),
+ a_delta(d_delta.begin(), p, false, true) {
}
predMod::predMod(Rcpp::NumericVector coef0)
- : d_coef0(coef0), d_delta(coef0.size()),
+ : d_coef(coef0.size()), d_coef0(coef0), d_delta(coef0.size()),
+ a_coef(d_coef.begin(), d_coef.size(), false, true),
a_coef0(d_coef0.begin(), d_coef0.size(), false, true),
a_delta(d_delta.begin(), d_delta.size(), false, true) {
}
+ void predMod::setCoef0(const Rcpp::NumericVector& cc)
+ throw (std::runtime_error) {
+ if (cc.size() != d_coef0.size())
+ throw std::runtime_error("size mismatch in setCoef0");
+ std::copy(cc.begin(), cc.end(), d_coef0.begin());
+ }
+
+ void predMod::installCoef0() {
+ std::copy(d_coef.begin(), d_coef.end(), d_coef0.begin());
+ }
+
densePred::densePred(Rcpp::NumericMatrix x)
throw (std::runtime_error)
: predMod(x.ncol()), d_X(x),
@@ -166,6 +196,18 @@
if (d_coef0.size() != d_X.ncol())
throw std::runtime_error("length(coef0) != ncol(X)");
}
+
+ void densePred::installCoef0() {
+ std::copy(d_coef.begin(), d_coef.end(), d_coef0.begin());
+ }
+
+ void densePred::setCoef0(const Rcpp::NumericVector& cc)
+ throw (std::runtime_error) {
+ if (cc.size() != d_coef0.size())
+ throw std::runtime_error("size mismatch in setCoef0");
+ std::copy(cc.begin(), cc.end(), d_coef0.begin());
+ }
+
// Returns the linear predictor from the full step
Rcpp::NumericVector
@@ -178,8 +220,10 @@
std::runtime_error("length(xwts) or length(wtres) != nrow(X)");
arma::vec a_xwts = arma::vec(xwts.begin(), n, false, true),
a_wtres = arma::vec(wtres.begin(), n, false, true);
- a_delta = solve(diagmat(a_xwts) * a_X, a_wtres);
- return NumericVector(wrap(a_X * (a_coef0 + a_delta)));
+ arma::vec tmp = solve(diagmat(a_xwts) * a_X, a_wtres);
+ std::copy(tmp.begin(), tmp.end(), a_delta.begin());
+ d_coef = d_coef0 + d_delta;
+ return NumericVector(wrap(a_X * a_coef));
}
// Returns the linear predictor from a step of fac
Modified: pkg/RcppModels/src/RcppModels.h
===================================================================
--- pkg/RcppModels/src/RcppModels.h 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/src/RcppModels.h 2010-12-01 02:08:30 UTC (rev 2635)
@@ -3,19 +3,19 @@
#include <Rcpp.h>
#include "glmFamily.h"
-class Irwls{
-public:
- const Rcpp::NumericMatrix x;
- const Rcpp::NumericVector y;
- arma::mat X;
- arma::vec Y;
- Rcpp::NumericVector wrt;
+// class Irwls{
+// public:
+// const Rcpp::NumericMatrix x;
+// const Rcpp::NumericVector y;
+// arma::mat X;
+// arma::vec Y;
+// Rcpp::NumericVector wrt;
- Irwls() {}
+// Irwls() {}
- Irwls(Rcpp::NumericMatrix, Rcpp::NumericVector);
- arma::vec fit(Rcpp::NumericVector);
-};
+// Irwls(Rcpp::NumericMatrix, Rcpp::NumericVector);
+// arma::vec fit(Rcpp::NumericVector);
+// };
namespace glm {
class lmResp {
@@ -49,7 +49,7 @@
class glmResp : public lmResp {
public:
- glmResp() {}
+// glmResp() {}
glmResp(Rcpp::List, Rcpp::NumericVector y)
throw (std::runtime_error);
glmResp(Rcpp::List, Rcpp::NumericVector y,
@@ -75,14 +75,19 @@
const Rcpp::NumericVector& weights() const {return d_weights;}
const Rcpp::NumericVector& wtres() const {return d_wtres;}
const Rcpp::NumericVector& y() const {return d_y;}
+ const std::string& family() const {return d_fam.fam();}
+ const std::string& link() const {return d_fam.lnk();}
// double deviance() const;
double residDeviance() const;
double updateMu(const Rcpp::NumericVector&);
double updateWts();
double wrss() const {return d_wrss;}
Rcpp::NumericVector devResid() const;
+ Rcpp::NumericVector muEta() const;
+ Rcpp::NumericVector sqrtWrkWt() const;
+ Rcpp::NumericVector variance() const;
Rcpp::NumericVector wrkResids() const;
-
+ Rcpp::NumericVector wrkResp() const;
protected:
glmFamily d_fam;
const Rcpp::NumericVector d_n;
@@ -91,15 +96,18 @@
class predMod {
protected:
- Rcpp::NumericVector d_coef0, d_delta;
- arma::vec a_coef0, a_delta;
+ Rcpp::NumericVector d_coef, d_coef0, d_delta;
+ arma::vec a_coef, a_coef0, a_delta;
public:
- predMod() {}
predMod(int);
predMod(Rcpp::NumericVector);
+
const Rcpp::NumericVector& delta() const {return d_delta;}
const Rcpp::NumericVector& coef0() const {return d_coef0;}
- void setCoef0(Rcpp::NumericVector);
+
+ void installCoef0();
+ void setCoef0(const Rcpp::NumericVector&)
+ throw (std::runtime_error);
};
class densePred : predMod {
@@ -107,15 +115,20 @@
Rcpp::NumericMatrix d_X;
arma::mat a_X;
public:
- densePred() {}
densePred(Rcpp::NumericMatrix) throw (std::runtime_error);
densePred(Rcpp::NumericMatrix, Rcpp::NumericVector)
throw (std::runtime_error);
const Rcpp::NumericMatrix& X() const {return d_X;}
- const arma::mat& aX() const {return a_X;}
+ const arma::mat& aX() const {return a_X;}
+ const Rcpp::NumericVector& coef() const {return d_coef;}
+ const Rcpp::NumericVector& coef0() const {return d_coef0;}
const Rcpp::NumericVector& delta() const {return d_delta;}
- const Rcpp::NumericVector& coef0() const {return d_coef0;}
+ void installCoef0();
+ void setCoef0(const Rcpp::NumericVector&)
+ throw (std::runtime_error);
+// Choose better names than gamma and dgamma. These conflict with
+// gamma functions, gamma distributions, the Gamma family, ...
Rcpp::NumericVector gamma(const Rcpp::NumericVector&,
const Rcpp::NumericVector&)
throw (std::runtime_error);
Modified: pkg/RcppModels/src/glmFamily.cpp
===================================================================
--- pkg/RcppModels/src/glmFamily.cpp 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/src/glmFamily.cpp 2010-12-01 02:08:30 UTC (rev 2635)
@@ -62,25 +62,25 @@
varFuncs["poisson"] = &identf; // x
}
- glmFamily::glmFamily() {
- if (!lnks.count("identity")) initMaps();
- }
+ // glmFamily::glmFamily() {
+ // if (!lnks.count("identity")) initMaps();
+ // }
- glmFamily::glmFamily(List ll) : lst(ll) {
+ glmFamily::glmFamily(List ll) throw (std::runtime_error)
+ : lst(ll) {
if (as<string>(lst.attr("class")) != "family")
- Rf_error("glmFamily only from list of (S3) class \"family\"");
-
- CharacterVector fam = lst["family"], llink = lst["link"];
- char *pt = fam[0]; family = string(pt);
- pt = llink[0]; link = string(pt);
+ throw std::runtime_error("glmFamily requires a list of (S3) class \"family\"");
+ CharacterVector ff = lst["family"], lnk = lst["link"];
+ d_family = as<std::string>(ff);
+ d_link = as<std::string>(lnk);
if (!lnks.count("identity")) initMaps();
}
Rcpp::NumericVector
glmFamily::linkFun(Rcpp::NumericVector const &mu) const {
- if (lnks.count(link)) { // sapply the known scalar function
- return NumericVector::import_transform(mu.begin(), mu.end(), lnks[link]);
+ if (lnks.count(d_link)) { // sapply the known scalar function
+ return NumericVector::import_transform(mu.begin(), mu.end(), lnks[d_link]);
} else { // use the R function
Function linkfun = ((const_cast<glmFamily*>(this))->lst)["linkfun"];
// The const_cast is needed so that this member function
@@ -92,8 +92,8 @@
Rcpp::NumericVector
glmFamily::linkInv(Rcpp::NumericVector const &eta) const {
- if (linvs.count(link)) {
- return NumericVector::import_transform(eta.begin(), eta.end(), linvs[link]);
+ if (linvs.count(d_link)) {
+ return NumericVector::import_transform(eta.begin(), eta.end(), linvs[d_link]);
} else {
Function linkinv = ((const_cast<glmFamily*>(this))->lst)["linkinv"];
return linkinv(eta);
@@ -102,8 +102,8 @@
Rcpp::NumericVector
glmFamily::muEta(Rcpp::NumericVector const &eta) const {
- if (muEtas.count(link)) {
- return NumericVector::import_transform(eta.begin(), eta.end(), muEtas[link]);
+ if (muEtas.count(d_link)) {
+ return NumericVector::import_transform(eta.begin(), eta.end(), muEtas[d_link]);
}
Function mu_eta = ((const_cast<glmFamily*>(this))->lst)["mu.eta"];
return mu_eta(eta);
@@ -111,8 +111,8 @@
Rcpp::NumericVector
glmFamily::variance(Rcpp::NumericVector const &mu) const {
- if (varFuncs.count(link)) {
- return NumericVector::import_transform(mu.begin(), mu.end(), varFuncs[link]);
+ if (varFuncs.count(d_link)) {
+ return NumericVector::import_transform(mu.begin(), mu.end(), varFuncs[d_link]);
}
Function vv = ((const_cast<glmFamily*>(this))->lst)["variance"];
return vv(mu);
@@ -122,8 +122,8 @@
glmFamily::devResid(Rcpp::NumericVector const &mu,
Rcpp::NumericVector const &weights,
Rcpp::NumericVector const &y) const {
- if (devRes.count(family)) {
- double (*f)(double, double, double) = devRes[family];
+ if (devRes.count(d_family)) {
+ double (*f)(double, double, double) = devRes[d_family];
int nobs = mu.size();
NumericVector ans(nobs);
double *mm = mu.begin(), *ww = weights.begin(),
Modified: pkg/RcppModels/src/glmFamily.h
===================================================================
--- pkg/RcppModels/src/glmFamily.h 2010-11-30 22:07:33 UTC (rev 2634)
+++ pkg/RcppModels/src/glmFamily.h 2010-12-01 02:08:30 UTC (rev 2635)
@@ -12,14 +12,17 @@
double(*)(double,double,double)> drmap;
class glmFamily {
- std::string family, link; // as in the R glm family
+ std::string d_family, d_link; // as in the R glm family
Rcpp::List lst; // original list from R
public:
- glmFamily();
+// glmFamily();
- glmFamily(Rcpp::List);
+ glmFamily(Rcpp::List) throw (std::runtime_error);
+ const std::string& fam() const {return d_family;}
+ const std::string& lnk() const {return d_link;}
+
// initialize the associative arrays of scalar functions
void initMaps();
Added: pkg/RcppModels/tests/doRUnit.R
===================================================================
--- pkg/RcppModels/tests/doRUnit.R (rev 0)
+++ pkg/RcppModels/tests/doRUnit.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,32 @@
+#### doRUnit.R --- Run RUnit tests
+####------------------------------------------------------------------------
+
+### borrowed from package Rcpp
+
+### MM: Vastly changed: This should also be "runnable" for *installed*
+## package which has no ./tests/
+## ----> put the bulk of the code e.g. in ../inst/unitTests/runTests.R :
+
+if( identical( .Platform$OS.type, "windows" ) && identical( .Platform$r_arch, "x64" ) ){
+ print( "unit tests not run on windows 64 (workaround alert)" )
+} else {
+ if(require("RUnit", quietly = TRUE)) {
+ pkg <- "RcppModels"
+
+ require( pkg, character.only=TRUE)
+
+ path <- system.file("unitTests", package = pkg)
+
+ stopifnot(file.exists(path), file.info(path.expand(path))$isdir)
+
+ # without this, we get unit test failures
+ Sys.setenv( R_TESTS = "" )
+
+ Rcpp.unit.test.output.dir <- getwd()
+
+ source(file.path(path, "runTests.R"), echo = TRUE)
+
+ } else {
+ print( "package RUnit not available, cannot run unit tests" )
+ }
+}
Added: pkg/RcppModels/tests/pois.R
===================================================================
--- pkg/RcppModels/tests/pois.R (rev 0)
+++ pkg/RcppModels/tests/pois.R 2010-12-01 02:08:30 UTC (rev 2635)
@@ -0,0 +1,22 @@
+library(RcppModels)
+
+## Dobson (1990) Page 93: Randomized Controlled Trial :
+d.AD <- data.frame(treatment = gl(3,3), outcome = gl(3,1,9),
+ counts = c(18,17,15,20,10,20,25,13,12))
+res <- fastGlm(counts ~ outcome + treatment, poisson, d.AD)
+pp <- res$pmod
+rr <- res$rmod
+rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres))
+(oldwrss <- rr$updateWts())
+pp$installCoef0()
+rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres))
+(oldwrss <- rr$updateWts())
+pp$installCoef0()
+(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
+(oldwrss <- rr$updateWts())
+pp$installCoef0()
+(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
+(oldwrss <- rr$updateWts())
+pp$installCoef0()
+(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
+
More information about the Rcpp-commits
mailing list