[Lme4-commits] r1502 - in pkg/lme4Eigen: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 10 21:50:42 CET 2012
Author: dmbates
Date: 2012-01-10 21:50:39 +0100 (Tue, 10 Jan 2012)
New Revision: 1502
Modified:
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
pkg/lme4Eigen/src/respModule.cpp
Log:
Got nlmer working, although it can take many iterations.
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-01-10 04:15:49 UTC (rev 1501)
+++ pkg/lme4Eigen/R/lmer.R 2012-01-10 20:50:39 UTC (rev 1502)
@@ -1,7 +1,5 @@
-##' Fit a linear mixed model or a generalized linear mixed model or a nonlinear mixed model.
+##' Fit a linear mixed-effects model
##'
-##'
-##'
##' @title Fit a linear mixed-effects model
##' @param formula a two-sided linear formula object describing the
##' fixed-effects part of the model, with the response on the left of a
@@ -84,7 +82,32 @@
if (devFunOnly) return(devfun)
if (verbose) control$iprint <- as.integer(verbose)
- opt <- bobyqa(reTrms$theta, devfun, reTrms$lower, control = control)
+
+ lower <- reTrms$lower
+ xst <- rep.int(0.1, length(lower))
+ nM <- NelderMead$new(lower=lower, upper=rep.int(Inf, length(lower)), xst=0.2*xst,
+ x0=rho$pp$theta, xt=xst*0.0001)
+ cc <- do.call(function(iprint=0L, maxfun=10000L, FtolAbs=1e-5,
+ FtolRel=1e-15, XtolRel=1e-7,
+ MinfMax=.Machine$double.xmin) {
+ if (length(list(...))>0) warning("unused control arguments ignored")
+ list(iprint=iprint, maxfun=maxfun, FtolAbs=FtolAbs, FtolRel=FtolRel,
+ XtolRel=XtolRel, MinfMax=MinfMax)
+ }, control)
+ nM$setMaxeval(cc$maxfun)
+ nM$setFtolAbs(cc$FtolAbs)
+ nM$setFtolRel(cc$FtolRel)
+ nM$setMinfMax(cc$MinfMax)
+ nM$setIprint(cc$iprint)
+ while ((nMres <- nM$newf(devfun(nM$xeval()))) == 0L) {}
+ if (nMres < 0L) {
+ if (nMres > -4L)
+ stop("convergence failure, code ", nMres, " in NelderMead")
+ else
+ warning("failure to converge in 1000 evaluations")
+ }
+ opt <- list(fval=nM$value(), pars=nM$xpos(), code=nMres)
+# opt <- bobyqa(reTrms$theta, devfun, reTrms$lower, control = control)
sqrLenU <- rho$pp$sqrL(1.)
wrss <- rho$resp$wrss()
pwrss <- wrss + sqrLenU
@@ -369,7 +392,7 @@
##' @param s Number of parameters in the nonlinear mean function (nlmer only)
##'
##' @return a list of Zt, Lambdat, Lind, theta, lower, flist and cnms
-mkReTrms <- function(bars, fr, s = 1L) {
+mkReTrms <- function(bars, fr) {
if (!length(bars))
stop("No random effects terms specified in formula")
stopifnot(is.list(bars), all(sapply(bars, is.language)),
@@ -442,7 +465,7 @@
thoff[i]))
}))))
thet <- numeric(sum(nth))
- ll <- list(Zt=Zt, theta=thet, Lind=as.integer(Lambdat at x),
+ ll <- list(Zt=Matrix::drop0(Zt), theta=thet, Lind=as.integer(Lambdat at x),
Gp=unname(c(0L, cumsum(nb))))
## lower bounds on theta elements are 0 if on diagonal, else -Inf
ll$lower <- -Inf * (thet + 1)
@@ -691,7 +714,7 @@
##' @param control a list of control parameters passed to the optimizer.
nlmer <- function(formula, data, control = list(), start = NULL, verbose = 0L,
nAGQ = 1L, subset, weights, na.action, offset,
- contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-10,
+ contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-5,
optimizer=c("NelderMead","bobyqa"), ...)
{
mf <- mc <- match.call()
@@ -737,7 +760,7 @@
for (nm in pnames) # convert these variables in fr to indicators
frE[[nm]] <- as.numeric(rep(nm == pnames, each = n))
# random-effects module
- reTrms <- mkReTrms(findbars(formula[[3]]), frE, s = s)
+ reTrms <- mkReTrms(findbars(formula[[3]]), frE)
fe.form <- nlform
fe.form[[3]] <- formula[[3]]
@@ -747,7 +770,7 @@
p <- ncol(X)
if ((qrX <- qr(X))$rank < p)
stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
- rho <- as.environment(list(verbose=verbose, tolPwrss=tolPwrss))
+ rho <- as.environment(list(verbose=verbose, tolPwrss=0.001))
parent.env(rho) <- parent.frame()
rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")],
list(X=X, S=s, Xwts = rr$sqrtXwt,
@@ -759,8 +782,40 @@
rho$lower <- reTrms$lower
devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
if (devFunOnly && !nAGQ) return(devfun)
+## FIXME: change this to something sensible for NelderMead
+
+ devfun(rho$pp$theta) # initial coarse evaluation to get u0 and beta0
+ rho$u0 <- rho$pp$u0
+ rho$beta0 <- rho$pp$beta0
+ rho$tolPwrss <- tolPwrss
control$iprint <- min(verbose, 3L)
- opt <- bobyqa(rho$pp$theta, devfun, rho$lower, control=control)
+ lower <- reTrms$lower
+ xst <- rep.int(0.1, length(lower))
+ nM <- NelderMead$new(lower=lower, upper=rep.int(Inf, length(lower)), xst=rep.int(-0.1, length(lower)),
+ x0=rho$pp$theta, xt=rep.int(0.00001, length(lower)))
+ cc <- do.call(function(iprint=0L, maxfun=10000L, FtolAbs=1e-5,
+ FtolRel=1e-15, XtolRel=1e-7,
+ MinfMax=.Machine$double.xmin) {
+ if (length(list(...))>0) warning("unused control arguments ignored")
+ list(iprint=iprint, maxfun=maxfun, FtolAbs=FtolAbs, FtolRel=FtolRel,
+ XtolRel=XtolRel, MinfMax=MinfMax)
+ }, control)
+ nM$setMaxeval(cc$maxfun)
+ nM$setFtolAbs(cc$FtolAbs)
+ nM$setFtolRel(cc$FtolRel)
+ nM$setMinfMax(cc$MinfMax)
+ nM$setIprint(cc$iprint)
+ while ((nMres <- nM$newf(devfun(nM$xeval()))) == 0L) {}
+ if (nMres < 0L) {
+ if (nMres > -4L)
+ stop("convergence failure, code ", nMres, " in NelderMead")
+ else
+ warning("failure to converge in ", cc$maxfun, " evaluations")
+ }
+ opt <- list(fval=nM$value(), pars=nM$xpos(), code=nMres)
+ cat(sprintf("After first opt: deviance = %g\n", nM$value()))
+ cat("theta: ", opt$pars, "\n")
+ cat("beta: ", rho$pp$beta0, "\n")
if (nAGQ > 0L) {
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
rho$u0 <- rho$pp$u0
@@ -802,7 +857,7 @@
if (nMres > -4L)
stop("convergence failure, code ", nMres, " in NelderMead")
else
- warning("failure to converge in 1000 evaluations")
+ warning("failure to converge in ", cc$maxfun, " evaluations")
}
list(fval=nM$value(), pars=nM$xpos(), code=nMres)
})
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-01-10 04:15:49 UTC (rev 1501)
+++ pkg/lme4Eigen/src/external.cpp 2012-01-10 20:50:39 UTC (rev 1502)
@@ -309,17 +309,13 @@
END_RCPP;
}
- static inline double prss(lmResp *rp, merPredD *pp, double fac) {
- return rp->wrss() + (fac ? pp->sqrL(fac) : pp->u0().squaredNorm());
- }
-
void nstepFac(nlsResp *rp, merPredD *pp, int verb) {
- double prss0(prss(rp, pp, 0.));
+ double prss0(pwrss(rp, pp, 0.));
for (double fac = 1.; fac > 0.001; fac /= 2.) {
double prss1 = rp->updateMu(pp->linPred(fac)) + pp->sqrL(fac);
if (verb > 3)
- ::Rprintf("pwrss0=%10g, diff=%10g, fac=%6.4f\n",
+ ::Rprintf("prss0=%10g, diff=%10g, fac=%6.4f\n",
prss0, prss0 - prss1, fac);
if (prss1 < prss0) {
pp->installPars(fac);
@@ -329,7 +325,6 @@
throw runtime_error("step factor reduced below 0.001 without reducing pwrss");
}
-#define MAXITER 30
static void prssUpdate(nlsResp *rp, merPredD *pp, int verb, bool uOnly, double tol) {
bool cvgd(false);
for (int it=0; it < MAXITER; ++it) {
@@ -337,13 +332,16 @@
pp->updateXwts(rp->sqrtXwt());
pp->updateDecomp();
pp->updateRes(rp->wtres());
- if ((uOnly ? pp->solveU() : pp->solve())/pwrss(rp, pp, 0.) < tol) {
+ double ccrit((uOnly ? pp->solveU() : pp->solve())/pwrss(rp, pp, 0.));
+ if (verb > 3)
+ ::Rprintf("ccrit=%10g, tol=%10g\n", ccrit, tol);
+ if (ccrit < tol) {
cvgd = true;
break;
}
nstepFac(rp, pp, verb);
}
- if (!cvgd) throw runtime_error("pwrss failed to converge in 30 iterations");
+ if (!cvgd) throw runtime_error("prss failed to converge in 30 iterations");
}
SEXP nlmerLaplace(SEXP pp_, SEXP rp_, SEXP theta_, SEXP u0_, SEXP beta0_,
@@ -352,7 +350,7 @@
XPtr<nlsResp> rp(rp_);
XPtr<merPredD> pp(pp_);
-
+// int verb(::Rf_asInteger(verbose_));
pp->setTheta(as<MVec>(theta_));
pp->setU0(as<MVec>(u0_));
pp->setBeta0(as<MVec>(beta0_));
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2012-01-10 04:15:49 UTC (rev 1501)
+++ pkg/lme4Eigen/src/predModule.cpp 2012-01-10 20:50:39 UTC (rev 1502)
@@ -52,10 +52,8 @@
d_VtV.setZero().selfadjointView<Eigen::Upper>().rankUpdate(d_V.adjoint());
d_RX.compute(d_VtV); // ensure d_RX is initialized even in the 0-column X case
-// Rcpp::Rcout << "before setTheta" << std::endl;
setTheta(d_theta); // starting values into Lambda
d_L.cholmod().final_ll = 1; // force an LL' decomposition
-// Rcpp::Rcout << "before updateLamtUt" << std::endl;
updateLamtUt();
d_L.analyzePattern(d_LamtUt); // perform symbolic analysis
if (d_L.info() != Eigen::Success)
@@ -67,7 +65,7 @@
// sparse/sparse matrix multiplication pruning zeros. The
// Cholesky decomposition croaks if the structure of d_LamtUt changes.
MVec(d_LamtUt._valuePtr(), d_LamtUt.nonZeros()).setZero();
- for (Index j = 0; j < d_Ut.cols(); ++j) {
+ for (Index j = 0; j < d_Ut.outerSize(); ++j) {
for(MSpMatrixd::InnerIterator rhsIt(d_Ut, j); rhsIt; ++rhsIt) {
Scalar y(rhsIt.value());
Index k(rhsIt.index());
@@ -80,7 +78,6 @@
}
}
}
-// Rcpp::Rcout << "end of updateLamtUt" << std::endl;
}
VectorXd merPredD::b(const double& f) const {return d_Lambdat.adjoint() * u(f);}
@@ -156,17 +153,15 @@
}
void merPredD::updateXwts(const VectorXd& sqrtXwt) {
-// Rcpp::Rcout << "in updateXwts" << std::endl;
if (d_Xwts.size() != sqrtXwt.size())
throw invalid_argument("updateXwts: dimension mismatch");
std::copy(sqrtXwt.data(), sqrtXwt.data() + sqrtXwt.size(), d_Xwts.data());
-// Rcpp::Rcout << "after copy, d_Xwts = " << d_Xwts.adjoint() << std::endl;
- if (sqrtXwt.size() == d_V.rows()) {
+ if (sqrtXwt.size() == d_V.rows()) { // W is diagonal
d_V = sqrtXwt.asDiagonal() * d_X;
for (int j = 0; j < d_N; ++j)
- for (MSpMatrixd::InnerIterator Uit(d_Ut, j), Zit(d_Zt, j);
- Uit && Zit; ++Uit, ++Zit)
- Uit.valueRef() = Zit.value() * sqrtXwt.data()[j];
+ for (MSpMatrixd::InnerIterator Utj(d_Ut, j), Ztj(d_Zt, j);
+ Utj && Ztj; ++Utj, ++Ztj)
+ Utj.valueRef() = Ztj.value() * sqrtXwt.data()[j];
} else {
SpMatrixd W(d_V.rows(), sqrtXwt.size());
const double *pt = sqrtXwt.data();
@@ -176,35 +171,22 @@
W.insertBack(j % d_V.rows(), j) = *pt;
}
W.finalize();
- // Rcpp::Rcout << "in nlmer block: W(" << W.rows() << "," << W.cols()
- // << "), outer: " << MiVec(W._outerIndexPtr(), W.outerSize()).adjoint()
- // << std::endl;
- // Rcpp::Rcout << "inner: " << MiVec(W._innerIndexPtr(), W.nonZeros()).adjoint()
- // << std::endl;
- // Rcpp::Rcout << "values: " << MVec(W._valuePtr(), W.nonZeros()).adjoint()
- // << std::endl;
d_V = W * d_X;
SpMatrixd Ut(d_Zt * W.adjoint());
- // Rcpp::Rcout << "Ut(" << Ut.rows() << "," << Ut.cols()
- // << "), outer: " << MiVec(Ut._outerIndexPtr(), Ut.outerSize()).adjoint()
- // << std::endl;
- // Rcpp::Rcout << "inner: " << MiVec(Ut._innerIndexPtr(), Ut.nonZeros()).adjoint()
- // << std::endl;
- // Rcpp::Rcout << "values: " << MVec(Ut._valuePtr(), Ut.nonZeros()).adjoint()
- // << std::endl;
if (Ut.cols() != d_Ut.cols())
throw std::runtime_error("Size mismatch in updateXwts");
// More complex code to handle the pruning of zeros
MVec(d_Ut._valuePtr(), d_Ut.nonZeros()).setZero();
- for (int j = 0; j < d_Ut.cols(); ++j) {
+ for (int j = 0; j < d_Ut.outerSize(); ++j) {
MSpMatrixd::InnerIterator lhsIt(d_Ut, j);
- SpMatrixd::InnerIterator rhsIt(Ut, j);
- Index k(rhsIt.index());
- while (lhsIt && lhsIt.index() != k) ++lhsIt;
- if (lhsIt.index() != k)
- throw std::runtime_error("Pattern mismatch in updateXwts");
- lhsIt.valueRef() = rhsIt.value();
+ for (SpMatrixd::InnerIterator rhsIt(Ut, j); rhsIt; ++rhsIt, ++lhsIt) {
+ Index k(rhsIt.index());
+ while (lhsIt && lhsIt.index() != k) ++lhsIt;
+ if (lhsIt.index() != k)
+ throw std::runtime_error("Pattern mismatch in updateXwts");
+ lhsIt.valueRef() = rhsIt.value();
+ }
}
}
d_VtV.setZero().selfadjointView<Eigen::Upper>().rankUpdate(d_V.adjoint());
Modified: pkg/lme4Eigen/src/respModule.cpp
===================================================================
--- pkg/lme4Eigen/src/respModule.cpp 2012-01-10 04:15:49 UTC (rev 1501)
+++ pkg/lme4Eigen/src/respModule.cpp 2012-01-10 20:50:39 UTC (rev 1502)
@@ -15,7 +15,6 @@
using Rcpp::as;
using std::copy;
- using std::string;
using std::invalid_argument;
typedef Eigen::Map<VectorXd> MVec;
@@ -153,7 +152,7 @@
double nlsResp::Laplace(double ldL2, double ldRX2, double sqrL) const {
double lnum = 2.* PI * (d_wrss + sqrL), n = d_y.size();
- return ldL2 + n * (1 + log(lnum / n));
+ return ldL2 + n * (1 + std::log(lnum / n));
}
double nlsResp::updateMu(const VectorXd& gamma) {
@@ -165,16 +164,16 @@
const double *gg = lp.data();
for (int p = 0; p < d_pnames.size(); ++p) {
- string pn(d_pnames[p]);
+ std::string pn(d_pnames[p]);
NumericVector pp = d_nlenv.get(pn);
- copy(gg + n * p, gg + n * (p + 1), pp.begin());
+ std::copy(gg + n * p, gg + n * (p + 1), pp.begin());
}
NumericVector rr = d_nlmod.eval(SEXP(d_nlenv));
if (rr.size() != n)
throw invalid_argument("dimension mismatch");
- copy(rr.begin(), rr.end(), d_mu.data());
+ std::copy(rr.begin(), rr.end(), d_mu.data());
NumericMatrix gr = rr.attr("gradient");
- copy(gr.begin(), gr.end(), d_sqrtXwt.data());
+ std::copy(gr.begin(), gr.end(), d_sqrtXwt.data());
return updateWrss();
}
More information about the Lme4-commits
mailing list