[Lme4-commits] r1572 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 8 00:18:50 CET 2012
Author: dmbates
Date: 2012-02-08 00:18:50 +0100 (Wed, 08 Feb 2012)
New Revision: 1572
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/glmFamily.cpp
pkg/lme4Eigen/src/glmFamily.h
Log:
Added an aic method to the glmFamily class but not yet working properly.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-02-07 19:43:59 UTC (rev 1571)
+++ pkg/lme4Eigen/R/AllClass.R 2012-02-07 23:18:50 UTC (rev 1572)
@@ -601,6 +601,16 @@
fields=list(Ptr="externalptr", family="family"),
methods=
list(
+ aic = function(y, n, mu, wt, dev) {
+ 'returns the value from the aic member function, which is actually the deviance'
+ nn <- length(y <- as.numeric(y))
+ stopifnot(length(n <- as.numeric(n)) == nn,
+ length(mu <- as.numeric(mu)) == nn,
+ length(wt <- as.numeric(wt)) == nn,
+ all(wt >= 0),
+ length(dev <- as.numeric(dev)) == 1L)
+ .Call(glmFamily_aic, ptr(), y, n, wt, mu, dev)
+ },
devResid = function(mu, weights, y) {
'applies the devResid function to mu, weights and y'
mu <- as.numeric(mu)
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-02-07 19:43:59 UTC (rev 1571)
+++ pkg/lme4Eigen/src/external.cpp 2012-02-07 23:18:50 UTC (rev 1572)
@@ -187,6 +187,16 @@
END_RCPP;
}
+ SEXP glmFamily_aic(SEXP ptr, SEXP y, SEXP n, SEXP mu, SEXP wt, SEXP dev) {
+ BEGIN_RCPP;
+ return ::Rf_ScalarReal(XPtr<glmFamily>(ptr)->aic(as<MVec>(y),
+ as<MVec>(n),
+ as<MVec>(mu),
+ as<MVec>(wt),
+ ::Rf_asReal(dev)));
+ END_RCPP;
+ }
+
SEXP glmFamily_muEta(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
return wrap(XPtr<glmFamily>(ptr)->muEta(as<MVec>(eta)));
@@ -816,22 +826,23 @@
CALLDEF(glm_wrkResids, 1),
CALLDEF(glm_wrkResp, 1),
- CALLDEF(glm_Laplace, 4), // methods
- CALLDEF(glm_updateMu, 2),
- CALLDEF(glm_updateWts, 1),
+ CALLDEF(glm_Laplace, 4), // methods
+ CALLDEF(glm_updateMu, 2),
+ CALLDEF(glm_updateWts, 1),
- CALLDEF(glmFamily_Create, 1), // generate external pointer
+ CALLDEF(glmFamily_Create, 1), // generate external pointer
- CALLDEF(glmFamily_link, 2), // methods
- CALLDEF(glmFamily_linkInv, 2),
+ CALLDEF(glmFamily_aic, 6), // methods
+ CALLDEF(glmFamily_link, 2),
+ CALLDEF(glmFamily_linkInv, 2),
CALLDEF(glmFamily_devResid, 4),
- CALLDEF(glmFamily_muEta, 2),
+ CALLDEF(glmFamily_muEta, 2),
CALLDEF(glmFamily_variance, 2),
- CALLDEF(glmerAGQ, 8),
- CALLDEF(glmerPwrssUpdate, 5),
- CALLDEF(glmerWrkIter, 2),
- CALLDEF(glmerLaplace, 8),
+ CALLDEF(glmerAGQ, 8),
+ CALLDEF(glmerPwrssUpdate, 5),
+ CALLDEF(glmerWrkIter, 2),
+ CALLDEF(glmerLaplace, 8),
CALLDEF(golden_Create, 2),
CALLDEF(golden_newf, 2),
Modified: pkg/lme4Eigen/src/glmFamily.cpp
===================================================================
--- pkg/lme4Eigen/src/glmFamily.cpp 2012-02-07 19:43:59 UTC (rev 1571)
+++ pkg/lme4Eigen/src/glmFamily.cpp 2012-02-07 23:18:50 UTC (rev 1572)
@@ -9,12 +9,13 @@
double glmFamily::epsilon = numeric_limits<double>::epsilon();
// initialize the function maps (i.e. associative arrays of functions)
- drmap glmFamily::devRes = drmap();
+ drmap glmFamily::devRes = drmap();
+ aicmap glmFamily::aics = aicmap();
- fmap glmFamily::linvs = fmap();
- fmap glmFamily::lnks = fmap();
- fmap glmFamily::muEtas = fmap();
- fmap glmFamily::varFuncs = fmap();
+ 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
@@ -69,7 +70,7 @@
// Functions from list components. This is a placeholder until I can
// do so.
d_devRes("c"), d_linkfun("c"), d_linkinv("c"),
- d_muEta("c"), d_variance("c") {
+ d_muEta("c"), d_variance("c"), d_aic("c") {
// d_devRes(wrap(ll["dev.resids"])),
// d_linkfun(wrap(ll["linkfun"])),
// d_linkinv(wrap(ll["linkinv"])),
@@ -78,14 +79,14 @@
if (!lst.inherits("family"))
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);
- d_linkinv = ll["linkinv"];
- d_linkfun = ll["linkfun"];
- d_muEta = ll["mu.eta"];
+ d_family = as<std::string>(ff);
+ d_link = as<std::string>(lnk);
+ d_linkinv = ll["linkinv"];
+ d_linkfun = ll["linkfun"];
+ d_muEta = ll["mu.eta"];
d_variance = ll["variance"];
- d_devRes = ll["dev.resids"];
-
+ d_devRes = ll["dev.resids"];
+ d_aic = ll["aic"];
if (!lnks.count("identity")) initMaps();
}
@@ -150,4 +151,16 @@
}
return ans;
}
+
+ double glmFamily::aic(const VectorXd& y, const VectorXd& n, const VectorXd& mu,
+ const VectorXd& wt, double dev) const {
+ int nn = mu.size();
+ if (aics.count(d_family)) return aics[d_family](y, n, mu, wt, dev);
+ SEXP ans = d_aic(NumericVector(y.data(), y.data() + nn),
+ NumericVector(n.data(), n.data() + nn),
+ NumericVector(mu.data(), mu.data() + nn),
+ NumericVector(wt.data(), wt.data() + nn),
+ ::Rf_ScalarReal(dev));
+ return Rcpp::as<double>(ans);
+ }
}
Modified: pkg/lme4Eigen/src/glmFamily.h
===================================================================
--- pkg/lme4Eigen/src/glmFamily.h 2012-02-07 19:43:59 UTC (rev 1571)
+++ pkg/lme4Eigen/src/glmFamily.h 2012-02-07 23:18:50 UTC (rev 1572)
@@ -4,6 +4,7 @@
#include <Rmath.h>
#include <RcppEigen.h>
+#include <cmath>
namespace glm {
using Eigen::VectorXd;
@@ -11,13 +12,15 @@
/** associative arrays of functions returning double from a double */
typedef std::map<std::string, double(*)(const double&)> fmap;
typedef std::map<std::string, double(*)(const double&,const double&,const double&)> drmap;
+ typedef std::map<std::string, double(*)(const VectorXd&,const VectorXd&,const VectorXd&,
+ const VectorXd&,double)> aicmap;
class glmFamily {
protected:
Rcpp::List lst; /**< original list from R */
std::string d_family, d_link; /**< as in the R glm family */
//@{ R functions from the family, as a fall-back
- Rcpp::Function d_devRes, d_linkfun, d_linkinv, d_muEta, d_variance;
+ Rcpp::Function d_devRes, d_linkfun, d_linkinv, d_muEta, d_variance, d_aic;
//@}
public:
glmFamily(Rcpp::List);
@@ -35,6 +38,10 @@
VectorXd devResid(const VectorXd&, const VectorXd&, const VectorXd&) const;
VectorXd muEta(const VectorXd&) const;
VectorXd variance(const VectorXd&) const;
+ double aic(const VectorXd&, const VectorXd&, const VectorXd&,
+ const VectorXd&, double) const;
+ /**< in keeping with the botched up nomenclature in the R glm function,
+ * the value of aic is the deviance */
private:
// Class members that are maps storing the scalar functions
static fmap
@@ -42,7 +49,8 @@
linvs, /**< scalar linkinv functions */
muEtas, /**< scalar muEta functions */
varFuncs; /**< scalar variance functions */
- static drmap devRes; /**< scalardeviance residuals functions */
+ static drmap devRes; /**< scalar deviance residuals functions */
+ static aicmap aics; /**< scalar aic functions */
static double epsilon; /**< Threshold for some comparisons */
@@ -126,6 +134,30 @@
return 2 * wt * (y_log_y(y, mu) - (y - mu));
}
//@}
+ //@{ AIC functions (which actually return the deviance, go figure)
+ struct CwiseRound {
+ const double operator()(const double& x) const {return nearbyint(x);}
+ };
+
+ static inline double
+ BinomialAIC (const VectorXd& y, const VectorXd& n, const VectorXd& mu,
+ const VectorXd& wt, double dev) {
+ Eigen::ArrayXd m((n.array() > 1).any() ? n : wt);
+ Eigen::ArrayXd yy((m * y.array()).unaryExpr(CwiseRound()));
+ m = m.unaryExpr(CwiseRound());
+ double ans(0.);
+ for (int i=0; i < mu.size(); ++i)
+ ans += (m[i] <= 0. ? 0. : (wt[i]/m[i]) * ::Rf_dbinom(yy[i], m[i], mu[i], true));
+ return (-2. * ans);
+ }
+ static inline double
+ PoissonAIC (const VectorXd& y, const VectorXd& n, const VectorXd& mu,
+ const VectorXd& wt, double dev) {
+ double ans(0.);
+ for (int i=0; i < mu.size(); ++i) ans += ::Rf_dpois(y[i], mu[i], true) * wt[i];
+ return (-2. * ans);
+ }
+ //@}
//@}
};
}
More information about the Lme4-commits
mailing list