[Rcpp-commits] r2470 - in pkg: . RcppModels RcppModels/man RcppModels/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 20 17:44:10 CET 2010
Author: dmbates
Date: 2010-11-20 17:44:09 +0100 (Sat, 20 Nov 2010)
New Revision: 2470
Added:
pkg/RcppModels/
pkg/RcppModels/DESCRIPTION
pkg/RcppModels/NAMESPACE
pkg/RcppModels/R/
pkg/RcppModels/man/
pkg/RcppModels/man/RccpModels-package.Rd
pkg/RcppModels/src/
pkg/RcppModels/src/Makevars
pkg/RcppModels/src/Makevars.win
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:
Initial version of linear, etc. models using RcppArmadillo
Added: pkg/RcppModels/DESCRIPTION
===================================================================
--- pkg/RcppModels/DESCRIPTION (rev 0)
+++ pkg/RcppModels/DESCRIPTION 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,13 @@
+Package: RcppModels
+Type: Package
+Title: Classes for linear and generalized linear and nonlinear models
+Version: 1.0
+Date: 2010-11-20
+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
+LinkingTo: RcppArmadillo, Rcpp
Added: pkg/RcppModels/NAMESPACE
===================================================================
--- pkg/RcppModels/NAMESPACE (rev 0)
+++ pkg/RcppModels/NAMESPACE 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,3 @@
+useDynLib(RcppModels)
+
+
Added: pkg/RcppModels/man/RccpModels-package.Rd
===================================================================
--- pkg/RcppModels/man/RccpModels-package.Rd (rev 0)
+++ pkg/RcppModels/man/RccpModels-package.Rd 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,27 @@
+\name{RcppModels-package}
+\alias{RcppModels-package}
+\alias{RcppModels}
+\docType{package}
+\title{
+ Reference classes for linear and generalized linear models
+}
+\description{
+ This package provides reference classes, defined through Rcpp Modules,
+ for statistical models defined by linear predictor expressions.
+}
+\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
+ }
+}
+\author{
+<Rcpp-core at r-forge.wu-wien.ac.at>
+}
+\keyword{ package }
+%\seealso{}
+%\examples{}
Added: pkg/RcppModels/src/Makevars
===================================================================
--- pkg/RcppModels/src/Makevars (rev 0)
+++ pkg/RcppModels/src/Makevars 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,5 @@
+## -*- mode: makefile; -*-
+## Use the R_HOME indirection to support installations of multiple R version
+
+PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+
Added: pkg/RcppModels/src/Makevars.win
===================================================================
--- pkg/RcppModels/src/Makevars.win (rev 0)
+++ pkg/RcppModels/src/Makevars.win 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,5 @@
+
+## This assume that we can call Rscript to ask Rcpp about its locations
+## Use the R_HOME indirection to support installations of multiple R version
+PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+
Added: pkg/RcppModels/src/Module.cpp
===================================================================
--- pkg/RcppModels/src/Module.cpp (rev 0)
+++ pkg/RcppModels/src/Module.cpp 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,103 @@
+#include "glmFamily.h"
+#include "RcppModels.h"
+
+using namespace Rcpp;
+
+RCPP_MODULE(RcppModels) {
+
+class_<Irwls>( "Irwls" )
+
+ .default_constructor()
+ .constructor(init_2<NumericMatrix,NumericVector>())
+
+ .field_readonly("x", &Irwls::x)
+ .field_readonly("y", &Irwls::y)
+
+ .method("fit", &Irwls::fit)
+ ;
+
+using namespace glm;
+
+class_<glmFamily>("glmFamily")
+
+ .default_constructor()
+ .constructor(init_1<List>())
+
+ .method("linkFun", &glmFamily::linkFun)
+ .method("linkInv", &glmFamily::linkInv)
+ .method("muEta", &glmFamily::muEta)
+ .method("devResid", &glmFamily::devResid)
+ .method("variance", &glmFamily::variance)
+ ;
+
+class_<lmResp>("lmResp")
+
+ .default_constructor()
+ .constructor(init_1<NumericVector>())
+ .constructor(init_2<NumericVector,NumericVector>())
+ .constructor(init_3<NumericVector,NumericVector,NumericVector>())
+
+ .property("mu", &lmResp::mu)
+ .property("offset", &lmResp::offset)
+ .property("sqrtXwt", &lmResp::sqrtXwt)
+ .property("sqrtrwt", &lmResp::sqrtrwt)
+ .property("weights", &lmResp::weights)
+ .property("wrss", &lmResp::wrss)
+ .property("wtres", &lmResp::wtres)
+ .property("y", &lmResp::y)
+
+ .method("updateMu", &lmResp::updateMu)
+ .method("updateWts", &lmResp::updateWts)
+ ;
+
+class_<glmResp>("glmResp")
+
+ .default_constructor()
+ .constructor(init_2<List,NumericVector>())
+ .constructor(init_3<List,NumericVector,NumericVector>())
+ .constructor(init_4<List,NumericVector,NumericVector,NumericVector>())
+ .constructor(init_5<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)
+
+ .method("devResid", &glmResp::devResid)
+ .method("residDeviance", &glmResp::residDeviance)
+ .method("updateMu", &glmResp::updateMu)
+ .method("updateWts", &glmResp::updateWts)
+ .method("wrkResids", &glmResp::wrkResids)
+ ;
+
+class_<predMod>("predMod")
+
+ .default_constructor()
+ .constructor(init_1<int>())
+ .constructor(init_1<NumericVector>())
+
+// For some reason this doesn't work.
+// .property("coef0", &predMod::coef0, &predMod::setCoef0)
+ .property("coef0", &predMod::coef0)
+ ;
+
+class_<densePred>("densePred")
+
+ .default_constructor()
+ .constructor(init_1<NumericMatrix>())
+ .constructor(init_2<NumericMatrix, NumericVector>())
+
+ .property("coef0", &densePred::coef0)
+ .property("delta", &densePred::delta)
+ .property("X", &densePred::X)
+
+ .method("gamma", &densePred::gamma)
+ ;
+
+}
Added: pkg/RcppModels/src/RcppModels.cpp
===================================================================
--- pkg/RcppModels/src/RcppModels.cpp (rev 0)
+++ pkg/RcppModels/src/RcppModels.cpp 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,189 @@
+#include "RcppModels.h"
+
+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);
+}
+
+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)
+ throw (std::runtime_error)
+ : d_y(y), d_weights(y.size(), 1.0), d_offset(y.size()),
+ d_mu(y.size()), d_sqrtrwt(y.size()),
+ d_sqrtXwt(y.size(), 1.0) {
+ init();
+ }
+
+ lmResp::lmResp(Rcpp::NumericVector y, Rcpp::NumericVector weights)
+ throw (std::runtime_error)
+ : d_y(y), d_weights(weights), d_offset(y.size()),
+ d_mu(y.size()), d_sqrtrwt(y.size()),
+ d_sqrtXwt(y.size(), 1) {
+ if (weights.size() != y.size())
+ throw std::runtime_error(
+ "lengths of y and wts must agree");
+ init();
+ }
+
+ lmResp::lmResp(Rcpp::NumericVector y, Rcpp::NumericVector weights,
+ Rcpp::NumericVector offset) throw (std::runtime_error)
+ : d_y(y), d_weights(weights), d_offset(offset),
+ d_mu(y.size()), d_sqrtrwt(y.size()),
+ d_sqrtXwt(y.size(), 1) {
+ int nn = y.size();
+ if (weights.size() != nn || offset.size() != nn)
+ throw std::runtime_error(
+ "lengths of y, weights and offset must agree");
+ init();
+ }
+
+ void lmResp::init() {
+ d_sqrtrwt = sqrt(d_weights);
+ std::copy(d_sqrtrwt.begin(), d_sqrtrwt.end(), d_sqrtXwt.begin());
+ }
+
+ double lmResp::updateMu(const Rcpp::NumericVector& gamma) {
+ d_mu = d_offset + gamma;
+ return updateWrss();
+ }
+
+ double lmResp::updateWrss() {
+ d_wtres = d_sqrtrwt * (d_y - d_mu);
+ d_wrss = std::inner_product(d_wtres.begin(), d_wtres.end(),
+ d_wtres.begin(), double());
+ return d_wrss;
+ }
+
+ glmResp::glmResp(Rcpp::List fam, Rcpp::NumericVector y)
+ throw (std::runtime_error)
+ : lmResp(y), d_fam(fam), d_eta(y.size()) {
+ }
+
+ glmResp::glmResp(Rcpp::List fam, Rcpp::NumericVector y,
+ Rcpp::NumericVector weights)
+ throw (std::runtime_error)
+ : lmResp(y, weights), d_fam(fam), d_eta(y.size()) {
+ }
+
+ glmResp::glmResp(Rcpp::List fam, Rcpp::NumericVector y,
+ Rcpp::NumericVector weights,
+ Rcpp::NumericVector offset)
+ throw (std::runtime_error)
+ : lmResp(y, weights, offset), d_fam(fam),
+ d_eta(y.size()) {
+ }
+
+ glmResp::glmResp(Rcpp::List fam, Rcpp::NumericVector y,
+ Rcpp::NumericVector weights,
+ Rcpp::NumericVector offset,
+ Rcpp::NumericVector n)
+ throw (std::runtime_error)
+ : lmResp(y, weights, offset), d_fam(fam), d_n(n),
+ d_eta(y.size()) {
+ if (n.size() != y.size())
+ throw std::runtime_error("lengths of y and n must agree");
+ }
+
+ glmResp::glmResp(Rcpp::List fam, Rcpp::NumericVector y,
+ Rcpp::NumericVector weights,
+ Rcpp::NumericVector offset,
+ Rcpp::NumericVector n, Rcpp::NumericVector eta)
+ throw (std::runtime_error)
+ : lmResp(y, weights, offset), d_fam(fam), d_n(n),
+ d_eta(eta) {
+ int nn = y.size();
+ if (n.size() != nn || eta.size() != nn )
+ throw std::runtime_error(
+ "lengths of y, n and eta must agree");
+ }
+
+ Rcpp::NumericVector glmResp::devResid() const {
+ return d_fam.devResid(d_mu, d_weights, d_y);
+ }
+
+ double glmResp::residDeviance() const {
+ return sum(devResid());
+ }
+
+ double glmResp::updateWts() {
+ d_sqrtrwt = sqrt(d_weights/d_fam.variance(d_mu));
+ NumericVector muEta = d_fam.muEta(d_eta);
+ std::transform(muEta.begin(), muEta.end(), d_sqrtrwt.begin(),
+ d_sqrtXwt.begin(), std::multiplies<double>());
+ return updateWrss();
+ }
+
+ Rcpp::NumericVector glmResp::wrkResids() const {
+ // This needs to be modified
+ return d_wtres;
+ }
+
+ double glmResp::updateMu(const Rcpp::NumericVector& gamma) {
+ d_eta = d_offset + gamma;
+ NumericVector mmu = d_fam.linkInv(d_eta);
+ std::copy(mmu.begin(), mmu.end(), d_mu.begin());
+ return updateWrss();
+ }
+
+ 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) {
+ }
+
+ predMod::predMod(Rcpp::NumericVector coef0)
+ : d_coef0(coef0), d_delta(coef0.size()),
+ a_coef0(d_coef0.begin(), d_coef0.size(), false, true),
+ a_delta(d_delta.begin(), d_delta.size(), false, true) {
+ }
+
+ densePred::densePred(Rcpp::NumericMatrix x)
+ throw (std::runtime_error)
+ : predMod(x.ncol()), d_X(x),
+ a_X(x.begin(), x.nrow(), x.ncol(), false, true) {
+ }
+
+ densePred::densePred(Rcpp::NumericMatrix x,
+ Rcpp::NumericVector coef0)
+ throw (std::runtime_error)
+ : predMod(coef0), d_X(x),
+ a_X(x.begin(), x.nrow(), x.ncol(), false, true) {
+ if (d_coef0.size() != d_X.ncol())
+ throw std::runtime_error("length(coef0) != ncol(X)");
+ }
+
+ // Returns the linear predictor from the full step
+ Rcpp::NumericVector
+ densePred::gamma(const Rcpp::NumericVector& xwts,
+ const Rcpp::NumericVector& wtres)
+ throw (std::runtime_error) {
+ int n = d_X.nrow();
+ if (xwts.size() != n || wtres.size() != n)
+ throw
+ 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)));
+ }
+
+ // Returns the linear predictor from a step of fac
+ Rcpp::NumericVector densePred::dgamma(double fac) {
+ return NumericVector(wrap(a_X * (a_coef0 + fac * a_delta)));
+ }
+}
Added: pkg/RcppModels/src/RcppModels.h
===================================================================
--- pkg/RcppModels/src/RcppModels.h (rev 0)
+++ pkg/RcppModels/src/RcppModels.h 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,124 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+#include <RcppArmadillo.h>
+#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;
+
+ Irwls() {}
+
+ Irwls(Rcpp::NumericMatrix, Rcpp::NumericVector);
+ arma::vec fit(Rcpp::NumericVector);
+};
+
+namespace glm {
+ class lmResp {
+ public:
+ lmResp() {}
+ lmResp(Rcpp::NumericVector y) throw (std::runtime_error);
+ lmResp(Rcpp::NumericVector y, Rcpp::NumericVector weights)
+ throw (std::runtime_error);
+ lmResp(Rcpp::NumericVector y, Rcpp::NumericVector weights,
+ Rcpp::NumericVector offset)
+ throw (std::runtime_error);
+
+ void init();
+ const Rcpp::NumericVector& mu() const {return d_mu;}
+ const Rcpp::NumericVector& offset() const {return d_offset;}
+ const Rcpp::NumericMatrix& sqrtXwt() const {return d_sqrtXwt;}
+ const Rcpp::NumericVector& sqrtrwt() const {return d_sqrtrwt;}
+ const Rcpp::NumericVector& weights() const {return d_weights;}
+ const Rcpp::NumericVector& wtres() const {return d_wtres;}
+ const Rcpp::NumericVector& y() const {return d_y;}
+ double wrss() const {return d_wrss;}
+ double updateMu(const Rcpp::NumericVector&);
+ double updateWts() {return updateWrss();}
+ double updateWrss();
+ protected:
+ double d_wrss;
+ const Rcpp::NumericVector d_y, d_weights, d_offset;
+ Rcpp::NumericVector d_mu, d_sqrtrwt, d_wtres;
+ Rcpp::NumericMatrix d_sqrtXwt;
+ };
+
+ class glmResp : public lmResp {
+ public:
+ glmResp() {}
+ glmResp(Rcpp::List, Rcpp::NumericVector y)
+ throw (std::runtime_error);
+ glmResp(Rcpp::List, Rcpp::NumericVector y,
+ Rcpp::NumericVector wts) throw (std::runtime_error);
+ glmResp(Rcpp::List, Rcpp::NumericVector y,
+ Rcpp::NumericVector wts,
+ Rcpp::NumericVector offset) throw (std::runtime_error);
+ glmResp(Rcpp::List, Rcpp::NumericVector y,
+ Rcpp::NumericVector wts,
+ Rcpp::NumericVector offset,
+ Rcpp::NumericVector n) throw (std::runtime_error);
+ glmResp(Rcpp::List, Rcpp::NumericVector y,
+ Rcpp::NumericVector wts,
+ Rcpp::NumericVector offset,
+ Rcpp::NumericVector n,
+ Rcpp::NumericVector eta) throw (std::runtime_error);
+
+ const Rcpp::NumericVector& eta() const {return d_eta;}
+ const Rcpp::NumericVector& mu() const {return d_mu;}
+ const Rcpp::NumericVector& offset() const {return d_offset;}
+ const Rcpp::NumericMatrix& sqrtXwt() const {return d_sqrtXwt;}
+ const Rcpp::NumericVector& sqrtrwt() const {return d_sqrtrwt;}
+ const Rcpp::NumericVector& weights() const {return d_weights;}
+ const Rcpp::NumericVector& wtres() const {return d_wtres;}
+ const Rcpp::NumericVector& y() const {return d_y;}
+// 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 wrkResids() const;
+
+ protected:
+ glmFamily d_fam;
+ const Rcpp::NumericVector d_n;
+ Rcpp::NumericVector d_eta;
+ };
+
+ class predMod {
+ protected:
+ Rcpp::NumericVector d_coef0, d_delta;
+ arma::vec 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);
+ };
+
+ class densePred : predMod {
+ protected:
+ 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 Rcpp::NumericVector& delta() const {return d_delta;}
+ const Rcpp::NumericVector& coef0() const {return d_coef0;}
+ Rcpp::NumericVector gamma(const Rcpp::NumericVector&,
+ const Rcpp::NumericVector&)
+ throw (std::runtime_error);
+ Rcpp::NumericVector dgamma(double fac);
+ };
+}
Added: pkg/RcppModels/src/glmFamily.cpp
===================================================================
--- pkg/RcppModels/src/glmFamily.cpp (rev 0)
+++ pkg/RcppModels/src/glmFamily.cpp 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,139 @@
+#include "glmFamily.h"
+#include <cmath>
+#include <limits>
+
+using namespace Rcpp;
+using namespace std;
+
+namespace glm {
+ // Establish the values for the class constants
+ double glmFamily::epsilon = numeric_limits<double>::epsilon();
+
+ // initialize the function maps (i.e. associative arrays of functions)
+ drmap glmFamily::devRes = drmap();
+
+ fmap glmFamily::linvs = fmap();
+ fmap glmFamily::lnks = fmap();
+ fmap glmFamily::muEtas = fmap();
+ fmap glmFamily::varFuncs = fmap();
+
+ void glmFamily::initMaps() {
+ // initialize the static maps. The identity link is
+ // guaranteed to be initialized if any maps are initialized
+ lnks["log"] = &log;
+ muEtas["log"] = linvs["log"] = &exp;
+
+ lnks["sqrt"] = &sqrt;
+ linvs["sqrt"] = &sqrf;
+ muEtas["sqrt"] = &twoxf;
+
+ lnks["identity"] = &identf;
+ linvs["identity"] = &identf;
+ muEtas["identity"] = &onef;
+
+ lnks["inverse"] = &inversef;
+ linvs["inverse"] = &inversef;
+ muEtas["inverse"] = &invderivf;
+
+ lnks["logit"] = &logitLink;
+ linvs["logit"] = &logitLinkInv;
+ muEtas["logit"] = &logitMuEta;
+
+ lnks["probit"] = &probitLink;
+ linvs["probit"] = &probitLinkInv;
+ muEtas["probit"] = &probitMuEta;
+
+// lnks["cloglog"] = &cloglogLink;
+ linvs["cloglog"] = &cloglogLinkInv;
+ muEtas["cloglog"] = &cloglogMuEta;
+
+ devRes["Gamma"] = &GammaDevRes;
+ varFuncs["Gamma"] = &sqrf; // x^2
+
+ devRes["binomial"] = &BinomialDevRes;
+ varFuncs["binomial"] = &x1mxf; // x * (1 - x)
+
+ devRes["gaussian"] = &GaussianDevRes;
+ varFuncs["gaussian"] = &onef; // 1
+
+ varFuncs["inverse.gaussian"] = &cubef; // x^3
+
+ devRes["poisson"] = &PoissonDevRes;
+ varFuncs["poisson"] = &identf; // x
+ }
+
+ glmFamily::glmFamily() {
+ if (!lnks.count("identity")) initMaps();
+ }
+
+ glmFamily::glmFamily(List ll) : 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);
+
+ 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]);
+ } else { // use the R function
+ Function linkfun = ((const_cast<glmFamily*>(this))->lst)["linkfun"];
+ // The const_cast is needed so that this member function
+ // can be const and also use the extraction of a list
+ // component.
+ return linkfun(mu);
+ }
+ }
+
+ Rcpp::NumericVector
+ glmFamily::linkInv(Rcpp::NumericVector const &eta) const {
+ if (linvs.count(link)) {
+ return NumericVector::import_transform(eta.begin(), eta.end(), linvs[link]);
+ } else {
+ Function linkinv = ((const_cast<glmFamily*>(this))->lst)["linkinv"];
+ return linkinv(eta);
+ }
+ }
+
+ Rcpp::NumericVector
+ glmFamily::muEta(Rcpp::NumericVector const &eta) const {
+ if (muEtas.count(link)) {
+ return NumericVector::import_transform(eta.begin(), eta.end(), muEtas[link]);
+ }
+ Function mu_eta = ((const_cast<glmFamily*>(this))->lst)["mu.eta"];
+ return mu_eta(eta);
+ }
+
+ Rcpp::NumericVector
+ glmFamily::variance(Rcpp::NumericVector const &mu) const {
+ if (varFuncs.count(link)) {
+ return NumericVector::import_transform(mu.begin(), mu.end(), varFuncs[link]);
+ }
+ Function vv = ((const_cast<glmFamily*>(this))->lst)["variance"];
+ return vv(mu);
+ }
+
+ Rcpp::NumericVector
+ 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];
+ int nobs = mu.size();
+ NumericVector ans(nobs);
+ double *mm = mu.begin(), *ww = weights.begin(),
+ *yy = y.begin(), *aa = ans.begin();
+ for (int i = 0; i < nobs; i++)
+ aa[i] = f(yy[i], mm[i], ww[i]);
+ return ans;
+ }
+ Function devres =
+ ((const_cast<glmFamily*>(this))->lst)["dev.resids"];
+ return devres(y, mu, weights);
+ }
+}
Added: pkg/RcppModels/src/glmFamily.h
===================================================================
--- pkg/RcppModels/src/glmFamily.h (rev 0)
+++ pkg/RcppModels/src/glmFamily.h 2010-11-20 16:44:09 UTC (rev 2470)
@@ -0,0 +1,127 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+#ifndef LME4_GLMFAMILY_H
+#define LME4_GLMFAMILY_H
+
+#include <Rcpp.h>
+#include <Rmath.h>
+
+namespace glm {
+ // Map (associative array) of functions returning double from a double
+ typedef std::map<std::string, double(*)(double)> fmap;
+ typedef std::map<std::string,
+ double(*)(double,double,double)> drmap;
+
+ class glmFamily {
+ std::string family, link; // as in the R glm family
+ Rcpp::List lst; // original list from R
+ public:
+
+ glmFamily();
+
+ glmFamily(Rcpp::List);
+
+ // initialize the associative arrays of scalar functions
+ void initMaps();
+
+ // Application of functions from the family
+ // The scalar transformations use compiled code when available
+ Rcpp::NumericVector linkFun(Rcpp::NumericVector const&) const;
+ Rcpp::NumericVector linkInv(Rcpp::NumericVector const&) const;
+ Rcpp::NumericVector devResid(
+ Rcpp::NumericVector const&,
+ Rcpp::NumericVector const&,
+ Rcpp::NumericVector const&) const;
+ Rcpp::NumericVector muEta(Rcpp::NumericVector const&) const;
+ Rcpp::NumericVector variance(Rcpp::NumericVector const&) const;
+ private:
+ // Class members that are maps storing the scalar functions
+ static fmap
+ lnks, // scalar link functions
+ linvs, // scalar linkinv functions
+ muEtas, // scalar muEta functions
+ varFuncs; // scalar variance functions
+ static drmap devRes;
+
+ // Thresholds common to the class
+ static double epsilon;
+
+ // Utility for deviance residuals
+ static inline double y_log_y(double y, double mu) {
+ return y ? y*log(y/mu) : 0;
+ }
+
+ // Scalar functions that will used in transform applications
+ static inline double cubef(double x) {return x * x * x;}
+ static inline double identf(double x) {return x;}
+ static inline double invderivf(double x) {return -1/(x * x);}
+ static inline double inversef(double x) {return 1/x;}
+ static inline double onef(double x) {return 1;}
+ static inline double sqrf(double x) {return x * x;}
+ static inline double twoxf(double x) {return 2 * x;}
+ static inline double x1mxf(double x) {
+ return std::max(epsilon, x * (1 - x));
+ }
+ static inline double logitLinkInv(double x) {
+ return Rf_plogis(x, 0., 1., 1, 0);
+ }
+ static inline double logitLink(double x) {
+ return Rf_qlogis(x, 0., 1., 1, 0);
+ }
+ static inline double logitMuEta(double x) {
+ return Rf_dlogis(x, 0., 1., 0);
+ }
+ static inline double probitLinkInv(double x) {
+ return Rf_pnorm5(x, 0., 1., 1, 0);
+ }
+ static inline double probitLink(double x) {
+ return Rf_qnorm5(x, 0., 1., 1, 0);
+ }
+ static inline double probitMuEta(double x) {
+ return Rf_dnorm4(x, 0., 1., 0);
+ }
+ // cumulative probability function of the complement of the
+ // Gumbel distribution (i.e. pgumbel(x,0.,1.,0) == 1-pgumbel2(-x,0.,1.,0))
+ static inline double
+ pgumbel2(double q, double loc, double scale, int lower_tail) {
+ q = (q - loc) / scale;
+ q = -exp(q);
+ return lower_tail ? -expm1(q) : exp(q);
+ }
+ static inline double cloglogLinkInv(double x) {
+ return pgumbel2(x, 0., 1., 1);
+ }
+
+ //density of the complement of the Gumbel distribution
+ static inline double
+ dgumbel2(double x, double loc, double scale, int give_log) {
+ x = (x - loc) / scale;
+ x = x - exp(x) - log(scale);
+ return give_log ? x : exp(x);
+ }
+ static inline double cloglogMuEta(double x) {
+ return dgumbel2(x, 0., 1., 0);
+ }
+
+ // deviance residuals functions
+ static inline double
+ BinomialDevRes(double y, double mu, double wt) {
+ return 2 * wt * (y_log_y(y, mu) + y_log_y(1 - y, 1 - mu));
+ }
+ static inline double
+ GammaDevRes (double y, double mu, double wt) {
+ return -2 * wt * y_log_y(y, mu) - (y - mu)/mu;
+ }
+ static inline double
+ GaussianDevRes(double y, double mu, double wt) {
+ double res = y - mu;
+ return wt * res * res;
+ }
+ static inline double
+ PoissonDevRes (double y, double mu, double wt) {
+ return 2 * wt * (y_log_y(y, mu) - (y - mu));
+ }
+ };
+}
+
+#endif /* LME4_GLMFAMILY_H */
+
More information about the Rcpp-commits
mailing list