[Lme4-commits] r1582 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 8 22:21:16 CET 2012
Author: dmbates
Date: 2012-02-08 22:21:15 +0100 (Wed, 08 Feb 2012)
New Revision: 1582
Modified:
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/glmFamily.cpp
pkg/lme4Eigen/src/glmFamily.h
pkg/lme4Eigen/src/respModule.cpp
pkg/lme4Eigen/src/respModule.h
Log:
Add the aic method for the glmFamily and glmResp classes. Still need to incorporate it in the evaluation of the deviance of the fitted model.
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/R/AllClass.R 2012-02-08 21:21:15 UTC (rev 1582)
@@ -462,6 +462,9 @@
n <<- if (!is.null(ll$n)) as.numeric(ll$n) else rep.int(1,length(y))
eta <<- numeric(length(y))
},
+ aic = function() {
+ .Call(glm_aic, ptr())
+ },
allInfo = function() {
'return all the information available on the object'
cbind(callSuper(),
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/src/external.cpp 2012-02-08 21:21:15 UTC (rev 1582)
@@ -78,6 +78,12 @@
END_RCPP;
}
+ SEXP glm_aic(SEXP ptr_) {
+ BEGIN_RCPP;
+ return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->aic());
+ END_RCPP;
+ }
+
SEXP glm_setN(SEXP ptr_, SEXP n) {
BEGIN_RCPP;
XPtr<glmResp>(ptr_)->setN(as<MVec>(n));
@@ -341,6 +347,7 @@
END_RCPP;
}
+
void nstepFac(nlsResp *rp, merPredD *pp, int verb) {
double prss0(pwrss(rp, pp, 0.));
@@ -816,15 +823,16 @@
CALLDEF(glm_setN, 2), // setters
- CALLDEF(glm_devResid, 1), // getters
- CALLDEF(glm_family, 1),
- CALLDEF(glm_link, 1),
- CALLDEF(glm_muEta, 1),
- CALLDEF(glm_resDev, 1),
- CALLDEF(glm_sqrtWrkWt, 1),
- CALLDEF(glm_variance, 1),
- CALLDEF(glm_wrkResids, 1),
- CALLDEF(glm_wrkResp, 1),
+ CALLDEF(glm_aic, 1), // getters
+ CALLDEF(glm_devResid, 1),
+ CALLDEF(glm_family, 1),
+ CALLDEF(glm_link, 1),
+ CALLDEF(glm_muEta, 1),
+ CALLDEF(glm_resDev, 1),
+ CALLDEF(glm_sqrtWrkWt, 1),
+ CALLDEF(glm_variance, 1),
+ CALLDEF(glm_wrkResids, 1),
+ CALLDEF(glm_wrkResp, 1),
CALLDEF(glm_Laplace, 4), // methods
CALLDEF(glm_updateMu, 2),
Modified: pkg/lme4Eigen/src/glmFamily.cpp
===================================================================
--- pkg/lme4Eigen/src/glmFamily.cpp 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/src/glmFamily.cpp 2012-02-08 21:21:15 UTC (rev 1582)
@@ -1,25 +1,146 @@
#include "glmFamily.h"
-#include <limits>
+#include <Rmath.h>
+#include <cmath>
using namespace Rcpp;
using namespace std;
namespace glm {
- // Establish the values for the class constants
- double glmFamily::epsilon = numeric_limits<double>::epsilon();
+ /**
+ * Utility to evaluate y * log(y/mu) with correct limit at y = 0
+ *
+ * @param y numerator of log argument, if non-zero
+ * @param mu denominator of log argument
+ *
+ * @return y * log(y/mu) with correct limit at at y = 0
+ */
+ static inline double y_log_y(const double& y, const double& mu) {
+ return y ? y * std::log(y/mu) : 0;
+ }
+
+ /** Cumulative probability function of the complement of the Gumbel distribution
+ *
+ * (i.e. pgumbel(q,0.,1.,0) == 1-pgumbel2(-q,0.,1.,0))
+ *
+ * @param q the quantile at which to evaluate the cumulative probability
+ * @param loc location parameter
+ * @param scale scale parameter
+ * @param lower_tail when zero evaluate the complement of the cdf
+ *
+ * @return Cumulative probability value or its complement, according to the value of lower_tail
+ */
+ static inline double
+ pgumbel2(const double& q, const double& loc, const double& scale, int lower_tail) {
+ double qq = (q - loc) / scale;
+ qq = -std::exp(qq);
+ return lower_tail ? -expm1(qq) : std::exp(qq);
+ }
+
+ /**
+ * density of the complement of the Gumbel distribution
+ *
+ * @param x numeric argument
+ * @param loc location parameter
+ * @param scale scale parameter
+ * @param give_log should the logarithm of the density be returned
+ *
+ * @return density or its logarithm, according to the value of give_log
+ */
+ static inline double
+ dgumbel2(const double& x, const double& loc, const double& scale, int give_log) {
+ double xx = (x - loc) / scale;
+ xx = xx - std::exp(xx) - std::log(scale);
+ return give_log ? xx : std::exp(xx);
+ }
+
+ template<typename Scalar>
+ struct Round {
+ const Scalar operator()(const Scalar& x) const {return nearbyint(x);}
+ };
+
+ //@{ Scalar functions for components
+ static inline double cubef(const double& x) {return x * x * x;}
+ static inline double expf(const double& x) {return std::exp(x);}
+ static inline double identf(const double& x) {return x;}
+ static inline double invderivf(const double& x) {return -1./(x * x);}
+ static inline double inversef(const double& x) {return 1./x;}
+ static inline double logf(const double& x) {return std::log(x);}
+ static inline double onef(const double& x) {return 1.;}
+ static inline double sqrf(const double& x) {return x * x;}
+ static inline double sqrtf(const double& x) {return std::sqrt(x);}
+ static inline double twoxf(const double& x) {return 2. * x;}
+ static inline double x1mxf(const double& x) {return std::max(epsilon, x*(1 - x));}
+ static inline double logitLinkInv(const double& x) {return ::Rf_plogis(x, 0., 1., 1, 0);}
+ static inline double logitLink(const double& x) {return ::Rf_qlogis(x, 0., 1., 1, 0);}
+ static inline double logitMuEta(const double& x) {return ::Rf_dlogis(x, 0., 1., 0);}
+ static inline double probitLinkInv(const double& x) {return ::Rf_pnorm5(x, 0., 1., 1, 0);}
+ static inline double probitLink(const double& x) {return ::Rf_qnorm5(x, 0., 1., 1, 0);}
+ static inline double probitMuEta(const double& x) {return ::Rf_dnorm4(x, 0., 1., 0);}
+ static inline double cloglogLinkInv(const double& x) {return pgumbel2(x, 0., 1., 1);}
+ //@}
- // initialize the function maps (i.e. associative arrays of functions)
- drmap glmFamily::devRes = drmap();
- aicmap glmFamily::aics = aicmap();
+ static inline double cloglogMuEta(const double& x) {return dgumbel2(x, 0., 1., 0);}
+
+ //@{ deviance residuals functions
+ static inline double
+ BinomialDevRes(const double& y, const double& mu, const double& wt) {
+ return 2 * wt * (y_log_y(y, mu) + y_log_y(1 - y, 1 - mu));
+ }
+ static inline double logYMu(const double& y, const double& mu) {
+ return y ? std::log(y/mu) : 0;
+ }
+ static inline double
+ GammaDevRes (const double& y, const double& mu, const double& wt) {
+ return -2 * wt * (logYMu(y, mu) - (y - mu)/mu);
+ }
+ static inline double
+ GaussianDevRes(const double& y, const double& mu, const double& wt) {
+ double res = y - mu;
+ return wt * res * res;
+ }
+ static inline double
+ PoissonDevRes (const double& y, const double& mu, const double& wt) {
+ return 2 * wt * (y_log_y(y, mu) - (y - mu));
+ }
+ //@}
+ //@{ AIC functions (which actually return the deviance, go figure)
+ 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(Round<double>()));
+ m = m.unaryExpr(Round<double>());
+ 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);
+ }
- fmap glmFamily::linvs = fmap();
- fmap glmFamily::lnks = fmap();
- fmap glmFamily::muEtas = fmap();
+ 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);
+ }
+ //@}
+
+ // initialize the function maps (i.e. associative arrays of
+ // functions). Needed because these are static maps.
+ drmap glmFamily::devRes = drmap();
+ aicmap glmFamily::aics = aicmap();
+
+ fmap glmFamily::linvs = fmap();
+ fmap glmFamily::lnks = fmap();
+ fmap glmFamily::muEtas = fmap();
fmap glmFamily::varFuncs = fmap();
+ /**
+ * Initialize the static maps. The identity link is guaranteed to be initialized if
+ * any of the maps are initialized. FIXME: Should probably check the size of the map instead?
+ *
+ */
void glmFamily::initMaps() {
- // initialize the static maps. The identity link is
- // guaranteed to be initialized if any maps are initialized
lnks["log"] = &logf;
muEtas["log"] = linvs["log"] = &expf;
@@ -50,6 +171,7 @@
devRes["Gamma"] = &GammaDevRes;
varFuncs["Gamma"] = &sqrf; // x^2
+ aics["binomial"] = &BinomialAIC;
devRes["binomial"] = &BinomialDevRes;
varFuncs["binomial"] = &x1mxf; // x * (1 - x)
@@ -58,6 +180,7 @@
varFuncs["inverse.gaussian"] = &cubef; // x^3
+ aics["poisson"] = &PoissonAIC;
devRes["poisson"] = &PoissonDevRes;
varFuncs["poisson"] = &identf; // x
}
Modified: pkg/lme4Eigen/src/glmFamily.h
===================================================================
--- pkg/lme4Eigen/src/glmFamily.h 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/src/glmFamily.h 2012-02-08 21:21:15 UTC (rev 1582)
@@ -2,10 +2,8 @@
#ifndef LME4_GLMFAMILY_H
#define LME4_GLMFAMILY_H
-#include <Rmath.h>
#include <RcppEigen.h>
-#include <cmath>
-
+#include <limits>
namespace glm {
using Eigen::VectorXd;
@@ -16,6 +14,12 @@
const VectorXd&,double)> aicmap;
class glmFamily {
+ public:
+
+ typedef std::map<std::string, double(*)(const double&)> fmap; /**< associative array of functions returning double from a double */
+ typedef std::map<std::string, double(*)(const double&,const double&,const double&)> drmap; /**< associative array of deviance residual functions */
+ typedef std::map<std::string, double(*)(const VectorXd&,const VectorXd&,const VectorXd&,
+ const VectorXd&,double)> aicmap;
protected:
Rcpp::List lst; /**< original list from R */
std::string d_family, d_link; /**< as in the R glm family */
@@ -52,114 +56,10 @@
static drmap devRes; /**< scalar deviance residuals functions */
static aicmap aics; /**< scalar aic functions */
- static double epsilon; /**< Threshold for some comparisons */
-
- /**
- * Utility to evaluate y * log(y/mu) with correct limit at y = 0
- *
- * @param y numerator of log argument, if non-zero
- * @param mu denominator of log argument
- *
- * @return y * log(y/mu) with correct limit at at y = 0
- */
- static inline double y_log_y(const double& y, const double& mu) {
- return y ? y * std::log(y/mu) : 0;
- }
-
-
- //@{ Scalar functions for components
- static inline double cubef(const double& x) {return x * x * x;}
- static inline double expf(const double& x) {return std::exp(x);}
- static inline double identf(const double& x) {return x;}
- static inline double invderivf(const double& x) {return -1./(x * x);}
- static inline double inversef(const double& x) {return 1./x;}
- static inline double logf(const double& x) {return std::log(x);}
- static inline double onef(const double& x) {return 1.;}
- static inline double sqrf(const double& x) {return x * x;}
- static inline double sqrtf(const double& x) {return std::sqrt(x);}
- static inline double twoxf(const double& x) {return 2. * x;}
- static inline double x1mxf(const double& x) {return std::max(epsilon, x*(1 - x));}
- static inline double logitLinkInv(const double& x) {return ::Rf_plogis(x, 0., 1., 1, 0);}
- static inline double logitLink(const double& x) {return ::Rf_qlogis(x, 0., 1., 1, 0);}
- static inline double logitMuEta(const double& x) {return ::Rf_dlogis(x, 0., 1., 0);}
- static inline double probitLinkInv(const double& x) {return ::Rf_pnorm5(x, 0., 1., 1, 0);}
- static inline double probitLink(const double& x) {return ::Rf_qnorm5(x, 0., 1., 1, 0);}
- static inline double probitMuEta(const 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(const double& q, const double& loc, const double& scale, int lower_tail) {
- double qq = (q - loc) / scale;
- qq = -std::exp(qq);
- return lower_tail ? -expm1(qq) : std::exp(qq);
- }
-
- static inline double cloglogLinkInv(const double& x) {
- return pgumbel2(x, 0., 1., 1);
- }
-
- //density of the complement of the Gumbel distribution
- static inline double
- dgumbel2(const double& x, const double& loc, const double& scale, int give_log) {
- double xx = (x - loc) / scale;
- xx = xx - std::exp(xx) - std::log(scale);
- return give_log ? xx : std::exp(xx);
- }
- static inline double cloglogMuEta(const double& x) {
- return dgumbel2(x, 0., 1., 0);
- }
-
- //@{ deviance residuals functions
- static inline double
- BinomialDevRes(const double& y, const double& mu, const double& wt) {
- return 2 * wt * (y_log_y(y, mu) + y_log_y(1 - y, 1 - mu));
- }
- static inline double logYMu(const double& y, const double& mu) {
- return y ? std::log(y/mu) : 0;
- }
- static inline double
- GammaDevRes (const double& y, const double& mu, const double& wt) {
- return -2 * wt * (logYMu(y, mu) - (y - mu)/mu);
- }
- static inline double
- GaussianDevRes(const double& y, const double& mu, const double& wt) {
- double res = y - mu;
- return wt * res * res;
- }
- static inline double
- PoissonDevRes (const double& y, const double& mu, const double& wt) {
- 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);
- }
- //@}
- //@}
};
+
+ static double epsilon(std::numeric_limits<double>::epsilon());
+ /**< Threshold for some comparisons */
}
#endif /* LME4_GLMFAMILY_H */
Modified: pkg/lme4Eigen/src/respModule.cpp
===================================================================
--- pkg/lme4Eigen/src/respModule.cpp 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/src/respModule.cpp 2012-02-08 21:21:15 UTC (rev 1582)
@@ -57,12 +57,12 @@
return d_wrss;
}
-/**
- * Set a new value of the offset.
- *
- * The values are copied into the d_offset member because it is mapped.
- * @param oo New value of the offset
- */
+ /**
+ * Set a new value of the offset.
+ *
+ * The values are copied into the d_offset member because that member is mapped.
+ * @param oo New value of the offset
+ */
void lmResp::setOffset(const VectorXd& oo) {
if (oo.size() != d_offset.size())
throw invalid_argument("setOffset: Size mismatch");
@@ -107,6 +107,10 @@
d_n(as<MVec>(n)) {
}
+ double glmResp::aic() const {
+ return d_fam.aic(d_y, d_n, d_mu, d_weights, resDev());
+ }
+
VectorXd glmResp::devResid() const {
return d_fam.devResid(d_mu, d_weights, d_y);
}
Modified: pkg/lme4Eigen/src/respModule.h
===================================================================
--- pkg/lme4Eigen/src/respModule.h 2012-02-08 20:27:07 UTC (rev 1581)
+++ pkg/lme4Eigen/src/respModule.h 2012-02-08 21:21:15 UTC (rev 1582)
@@ -102,6 +102,7 @@
const std::string& family() const {return d_fam.fam();}
const std::string& link() const {return d_fam.lnk();}
+ double aic() const;
double Laplace(double,double,double) const;
double resDev() const;
double updateMu(const Eigen::VectorXd&);
More information about the Lme4-commits
mailing list