[Lme4-commits] r1629 - in pkg/lme4Eigen: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 29 00:26:26 CET 2012
Author: dmbates
Date: 2012-02-29 00:26:26 +0100 (Wed, 29 Feb 2012)
New Revision: 1629
Added:
pkg/lme4Eigen/R/mcmcsamp.R
pkg/lme4Eigen/man/mcmcsamp.Rd
pkg/lme4Eigen/src/mcmcsamp.cpp
pkg/lme4Eigen/src/mcmcsamp.h
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/AllGeneric.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/predModule.cpp
pkg/lme4Eigen/src/predModule.h
Log:
Adding the shell of an mcmcsamp function.
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/DESCRIPTION 2012-02-28 23:26:26 UTC (rev 1629)
@@ -53,3 +53,4 @@
'predict.R'
'bootMer.R'
'nbinom.R'
+ 'mcmcsamp.R'
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/NAMESPACE 2012-02-28 23:26:26 UTC (rev 1629)
@@ -112,6 +112,7 @@
S3method(isREML,merMod)
S3method(logLik,merMod)
S3method(log,thpr)
+S3method(mcmcsamp,merMod)
S3method(model.frame,merMod)
S3method(model.matrix,merMod)
S3method(nobs,merMod)
Modified: pkg/lme4Eigen/R/AllGeneric.R
===================================================================
--- pkg/lme4Eigen/R/AllGeneric.R 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/R/AllGeneric.R 2012-02-28 23:26:26 UTC (rev 1629)
@@ -18,6 +18,19 @@
eval(mCall, parent.frame())
}
+##' Create a Markov chain Monte Carlo sample from the posterior
+##' distribution of the parameters
+##'
+##'
+##' @title Create an McMC sample
+##' @param object a fitted model object
+##' @param n number of samples to generate. Defaults to 1 but that
+##' value doesn't make sense.
+##' @param verbose should verbose output be given?
+##' @param ...
+##' @return a Markov chain Monte Carlo sample as a matrix
+mcmcsamp <- function(object, n = 1L, verbose = FALSE, ...) UseMethod("mcmcsamp")
+
##' Extract the residual standard error from a fitted model.
##'
##' This is a generic function. At present the only methods are for mixed-effects
Added: pkg/lme4Eigen/R/mcmcsamp.R
===================================================================
--- pkg/lme4Eigen/R/mcmcsamp.R (rev 0)
+++ pkg/lme4Eigen/R/mcmcsamp.R 2012-02-28 23:26:26 UTC (rev 1629)
@@ -0,0 +1,24 @@
+#' @S3method mcmcsamp merMod
+mcmcsamp.merMod <- function(object, n=1L, verbose=FALSE, saveb=FALSE, ...) {
+ n <- max(1L, as.integer(n)[1])
+ dd <- getME(object, "devcomp")$dims
+ ranef <- matrix(numeric(0), nrow = dd[["q"]], ncol = 0)
+ if (saveb) ranef <- matrix(, nrow = dd[["q"]], ncol = n)
+ sigma <- matrix(unname(sigma(object)), nrow = 1,
+ ncol = (if (dd[["useSc"]]) n else 0))
+ ff <- fixef(object)
+ fixef <- matrix(ff, nrow=dd[["p"]], ncol=n)
+ rownames(fixef) <- names(ff)
+ ## FIXME create a copy of the resp and pred modules
+ ans <- new("merMCMC",
+ Gp = object at Gp,
+ ST = matrix(.Call(mer_ST_getPars, object), dd[["np"]], n),
+ call = object at call,
+ dims = object at dims,
+ deviance = rep.int(unname(object at deviance[["ML"]]), n),
+ fixef = fixef,
+ nc = sapply(object at ST, nrow),
+ ranef = ranef,
+ sigma = sigma)
+ .Call(mer_MCMCsamp, ans, object)
+}
Added: pkg/lme4Eigen/man/mcmcsamp.Rd
===================================================================
--- pkg/lme4Eigen/man/mcmcsamp.Rd (rev 0)
+++ pkg/lme4Eigen/man/mcmcsamp.Rd 2012-02-28 23:26:26 UTC (rev 1629)
@@ -0,0 +1,24 @@
+\name{mcmcsamp}
+\alias{mcmcsamp}
+\title{Create an McMC sample}
+\usage{
+ mcmcsamp(object, n = 1L, verbose = FALSE, ...)
+}
+\arguments{
+ \item{object}{a fitted model object}
+
+ \item{n}{number of samples to generate. Defaults to 1
+ but that value doesn't make sense.}
+
+ \item{verbose}{should verbose output be given?}
+
+ \item{...}{}
+}
+\value{
+ a Markov chain Monte Carlo sample as a matrix
+}
+\description{
+ Create a Markov chain Monte Carlo sample from the
+ posterior distribution of the parameters
+}
+
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/src/external.cpp 2012-02-28 23:26:26 UTC (rev 1629)
@@ -276,7 +276,7 @@
}
static ArrayXd devcCol(const MiVec& fac, const ArrayXd& u, const ArrayXd& devRes) {
- ArrayXd ans(u * u);
+ ArrayXd ans(u.square());
for (int i = 0; i < devRes.size(); ++i) ans[fac[i] - 1] += devRes[i];
return ans;
}
@@ -294,7 +294,7 @@
throw std::invalid_argument("size of fac must match dimension of response vector");
pp->setTheta(as<MVec>(theta_));
- pp->setU0(as<MVec>(u0_));
+ pp->setU0( as<MVec>(u0_));
pp->setBeta0(as<MVec>(beta0_));
pwrssUpdate(rp, pp, 0, true, ::Rf_asReal(tol_)); // should be a no-op
const ArrayXd u0(pp->u0());
Added: pkg/lme4Eigen/src/mcmcsamp.cpp
===================================================================
--- pkg/lme4Eigen/src/mcmcsamp.cpp (rev 0)
+++ pkg/lme4Eigen/src/mcmcsamp.cpp 2012-02-28 23:26:26 UTC (rev 1629)
@@ -0,0 +1,258 @@
+//
+// mcmcsamp.cpp: implementation of mcmcsamp and related classes using Eigen
+//
+// Copyright (C) 2011-2012 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of lme4Eigen.
+
+#include "mcmcsamp.h"
+
+namespace lme4 {
+ using Rcpp::as;
+
+ static inline double pwrss(lme4Eigen::merPredD *pred, lme4Eigen::lmResp *resp) {
+ return pred->sqrL(1.) + resp->wrss();
+ }
+
+ static inline double sigmaML(lme4Eigen::merPredD *pred, lme4Eigen::lmResp *resp) {
+ return std::sqrt(pwrss(pred, resp)/double(resp->y().size()));
+ }
+
+ mcmcsamp::mcmcsamp(lme4Eigen::merPredD *pred, lme4Eigen::lmResp *resp,
+ SEXP dev, SEXP fixef, SEXP sigma, SEXP ranef)
+ : d_dev( as<MVec>(dev)),
+ d_fixef(as<MMat>(fixef)),
+ d_sigma(as<MVec>(sigma)),
+ d_ranef(as<MMat>(ranef))
+ {
+ Rcpp::RNGScope scope; // handles the getRNGstate/putRNGstate
+ bool sig( d_sigma.size() > 0);
+ bool rr( d_ranef.rows() > 0);
+ int n( resp->y().size());
+ int nsamp(d_dev.size());
+ int nth( pred->theta().size());
+ int p( pred->beta0().size());
+ int q( pred->u0().size());
+ double npq( n + q);
+ double lsigma(sig ? sigmaML(pred, resp) : 1.);
+ double wrss;
+
+ if (d_fixef.cols() != nsamp || d_fixef.rows() != p || (sig && d_sigma.size() != nsamp) ||
+ (ranef && (d_ranef.cols() != nsamp || d_ranef.rows() != p)))
+ throw std::invalid_argument("dimension mismatch");
+ if (nth > 1) ::Rf_error("only handling the simple (nth == 1) cases now");
+
+ for (int k = 0; k < nsamp; ++k) {
+ pred->updateDecomp();
+ pred->solve();
+ pred->MCMC_beta_u(lsigma);
+ d_fixef.col(k) = pred->beta(1.);
+ if (rr) d_ranef.col(k) = pred->b(1.);
+ wrss = resp->updateMu(pred->linPred(1.));
+ if (sig) d_sigma[k] = lsigma = std::sqrt((pred->sqrL(1.) + wrss)/::Rf_rchisq(npq));
+ }
+ }
+
+#if 0
+/**
+ * Generate a Markov-chain Monte Carlo sample from an mer object
+ *
+ * @param x pointer to an merMCMC object
+ * @param fm pointer to an mer object
+ *
+ * @return x with samples filled in
+ */
+ SEXP mer_MCMCsamp(SEXP , SEXP fm)
+ {
+ SEXP devsamp = GET_SLOT(x, lme4_devianceSym);
+ int *dims = DIMS_SLOT(x), nsamp = LENGTH(devsamp);
+ int n = dims[n_POS], np = dims[np_POS],
+ p = dims[p_POS], q = dims[q_POS];
+ double
+ *STsamp = REAL(GET_SLOT(x, lme4_STSym)),
+ *d = DEV_SLOT(fm), *dev = REAL(devsamp),
+ *sig = SLOT_REAL_NULL(x, lme4_sigmaSym),
+ *fixsamp = FIXEF_SLOT(x), *resamp = RANEF_SLOT(x);
+
+ GetRNGstate();
+ /* The first column of storage slots contains the fitted values */
+ for (int i = 1; i < nsamp; i++) {
+/* FIXME: This is probably wrong for a model with weights. */
+ if (sig) /* update and store sigma */
+ sig[i] = sqrt(d[pwrss_POS]/rchisq((double)(n + q)));
+ /* update L, RX, beta, u, eta, mu, res and d */
+ MCMC_beta_u(fm, sig ? sig[i] : 1, fixsamp + i * p,
+ resamp + (resamp ? i : 0) * q);
+ dev[i] = d[ML_POS];
+ /* update theta_T, theta_S and A */
+ MCMC_T(fm, sig ? sig[i] : 1);
+ MCMC_S(fm, sig ? sig[i] : 1);
+ ST_getPars(fm, STsamp + i * np); /* record theta */
+ }
+ PutRNGstate();
+ /* Restore pars from the first columns of the samples */
+ Memcpy(FIXEF_SLOT(fm), fixsamp, p);
+ ST_setPars(fm, STsamp);
+ update_ranef(fm);
+
+ return x;
+ }
+
+/**
+ * Update the fixed effects and the orthogonal random effects in an MCMC sample
+ * from an mer object.
+ *
+ * @param x an mer object
+ * @param sigma current standard deviation of the per-observation
+ * noise terms.
+ * @param fvals pointer to memory in which to store the updated beta
+ * @param rvals pointer to memory in which to store the updated b (may
+ * be (double*)NULL)
+ */
+ static void MCMC_beta_u(SEXP x, double sigma, double *fvals, double *rvals)
+ {
+ int *dims = DIMS_SLOT(x);
+ int i1 = 1, p = dims[p_POS], q = dims[q_POS];
+ double *V = V_SLOT(x), *fixef = FIXEF_SLOT(x), *muEta = MUETA_SLOT(x),
+ *u = U_SLOT(x), mone[] = {-1,0}, one[] = {1,0};
+ CHM_FR L = L_SLOT(x);
+ double *del1 = Calloc(q, double), *del2 = Alloca(p, double);
+ CHM_DN sol, rhs = N_AS_CHM_DN(del1, q, 1);
+ R_CheckStack();
+
+ if (V || muEta) {
+ error(_("Update not yet written"));
+ } else { /* Linear mixed model */
+ update_L(x);
+ update_RX(x);
+ lmm_update_fixef_u(x);
+ /* Update beta */
+ for (int j = 0; j < p; j++) del2[j] = sigma * norm_rand();
+ F77_CALL(dtrsv)("U", "N", "N", &p, RX_SLOT(x), &p, del2, &i1);
+ for (int j = 0; j < p; j++) fixef[j] += del2[j];
+ /* Update u */
+ for (int j = 0; j < q; j++) del1[j] = sigma * norm_rand();
+ F77_CALL(dgemv)("N", &q, &p, mone, RZX_SLOT(x), &q,
+ del2, &i1, one, del1, &i1);
+ sol = M_cholmod_solve(CHOLMOD_Lt, L, rhs, &c);
+ for (int j = 0; j < q; j++) u[j] += ((double*)(sol->x))[j];
+ M_cholmod_free_dense(&sol, &c);
+ update_mu(x); /* and parts of the deviance slot */
+ }
+ Memcpy(fvals, fixef, p);
+ if (rvals) {
+ update_ranef(x);
+ Memcpy(rvals, RANEF_SLOT(x), q);
+ }
+ Free(del1);
+ }
+
+/**
+ * Update the theta_T parameters from the ST arrays in place.
+ *
+ * @param x an mer object
+ * @param sigma current standard deviation of the per-observation
+ * noise terms.
+ */
+/* FIXME: Probably should fold this function into MCMC_S */
+ static void MCMC_T(SEXP x, double sigma)
+ {
+ int *Gp = Gp_SLOT(x), nt = (DIMS_SLOT(x))[nt_POS];
+ double **st = Alloca(nt, double*);
+ int *nc = Alloca(nt, int), *nlev = Alloca(nt, int);
+ R_CheckStack();
+
+ if (ST_nc_nlev(GET_SLOT(x, lme4_STSym), Gp, st, nc, nlev) < 2) return;
+ error("Code for non-trivial theta_T not yet written");
+ }
+
+/**
+ * Update the theta_S parameters from the ST arrays in place.
+ *
+ * @param x an mer object
+ * @param sigma current standard deviation of the per-observation
+ * noise terms.
+ */
+ static void MCMC_S(SEXP x, double sigma)
+ {
+ CHM_SP A = A_SLOT(x), Zt = Zt_SLOT(x);
+ int *Gp = Gp_SLOT(x), *ai = (int*)(A->i),
+ *ap = (int*)(A->p), *dims = DIMS_SLOT(x), *perm = PERM_VEC(x);
+ int annz = ap[A->ncol], info, i1 = 1, n = dims[n_POS],
+ nt = dims[nt_POS], ns, p = dims[p_POS], pos,
+ q = dims[q_POS], znnz = ((int*)(Zt->p))[Zt->ncol];
+ double *R, *ax = (double*)(A->x), *b = RANEF_SLOT(x),
+ *eta = ETA_SLOT(x), *offset = OFFSET_SLOT(x),
+ *rr, *ss, one = 1, *u = U_SLOT(x), *y = Y_SLOT(x);
+ int *nc = Alloca(nt, int), *nlev = Alloca(nt, int),
+ *spt = Alloca(nt + 1, int);
+ double **st = Alloca(nt, double*);
+ R_CheckStack();
+
+ ST_nc_nlev(GET_SLOT(x, lme4_STSym), Gp, st, nc, nlev);
+ ns = 0; /* ns is length(theta_S) */
+ spt[0] = 0; /* pointers into ss for terms */
+ for (int i = 0; i < nt; i++) {
+ ns += nc[i];
+ spt[i + 1] = spt[i] + nc[i];
+ }
+
+ if (annz == znnz) { /* Copy Z' to A unless A has new nonzeros */
+ Memcpy(ax, (double*)(Zt->x), znnz);
+ } else error("Code not yet written for MCMC_S with NLMMs");
+ /* Create T'Zt in A */
+ Tt_Zt(A, Gp, nc, nlev, st, nt);
+ /* Create P'u in ranef slot */
+ for (int i = 0; i < q; i++) b[perm[i]] = u[i];
+ /* Create X\beta + offset in eta slot */
+ for (int i = 0; i < n; i++) eta[i] = offset ? offset[i] : 0;
+ F77_CALL(dgemv)("N", &n, &p, &one, X_SLOT(x), &n,
+ FIXEF_SLOT(x), &i1, &one, eta, &i1);
+ /* Allocate R, rr and ss */
+ R = Alloca(ns * ns, double); /* crossproduct matrix then factor */
+ rr = Alloca(ns, double); /* row of model matrix for theta_S */
+ ss = Alloca(ns, double); /* right hand side, then theta_S */
+ R_CheckStack();
+ AZERO(R, ns * ns);
+ AZERO(ss, ns);
+ /* Accumulate crossproduct from pseudo-data part of model matrix */
+ for (int i = 0; i < q; i++) {
+ int sj = theta_S_ind(i, nt, Gp, nlev, spt);
+ AZERO(rr, ns);
+ rr[sj] = b[i];
+ F77_CALL(dsyr)("U", &ns, &one, rr, &i1, R, &ns);
+ }
+ /* Accumulate crossproduct and residual product of the model matrix. */
+ /* This is done one row at a time. Rows of the model matrix
+ * correspond to columns of T'Zt */
+ for (int j = 0; j < n; j++) { /* jth column of T'Zt */
+ AZERO(rr, ns);
+ for (int p = ap[j]; p < ap[j + 1]; p++) {
+ int i = ai[p]; /* row in T'Zt */
+ int sj = theta_S_ind(i, nt, Gp, nlev, spt);
+
+ rr[sj] += ax[p] * b[i];
+ ss[sj] += rr[sj] * (y[j] - eta[j]);
+ }
+ F77_CALL(dsyr)("U", &ns, &one, rr, &i1, R, &ns);
+ }
+ F77_CALL(dposv)("U", &ns, &i1, R, &ns, ss, &ns, &info);
+ if (info)
+ error(_("Model matrix for theta_S is not positive definite, %d."), info);
+ for (int j = 0; j < ns; j++) rr[j] = sigma * norm_rand();
+ /* Sample from the conditional Gaussian distribution */
+ F77_CALL(dtrsv)("U", "N", "N", &ns, R, &ns, rr, &i1);
+ for (int j = 0; j < ns; j++) ss[j] += rr[j];
+ /* Copy positive part of solution onto diagonals of ST */
+ pos = 0;
+ for (int i = 0; i < nt; i++) {
+ for (int j = 0; j < nc[i]; j++) {
+ st[i][j * (nc[i] + 1)] = (ss[pos] > 0) ? ss[pos] : 0;
+ pos++;
+ }
+ }
+ update_A(x);
+ }
+
+#endif
+}
Added: pkg/lme4Eigen/src/mcmcsamp.h
===================================================================
--- pkg/lme4Eigen/src/mcmcsamp.h (rev 0)
+++ pkg/lme4Eigen/src/mcmcsamp.h 2012-02-28 23:26:26 UTC (rev 1629)
@@ -0,0 +1,40 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+//
+// mcmcsamp.h: Markov-chain Monte Carlo sample class using Eigen
+//
+// Copyright (C) 2012 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of lme4.
+#ifndef LME4_MCMCSAMP_H
+#define LME4_MCMCSAMP_H
+
+#include "predModule.h"
+#include "respModule.h"
+
+namespace lme4 {
+ class mcmcsamp {
+ public:
+ typedef Eigen::ArrayXd Ar1;
+ typedef Eigen::Map<Ar1> MAr1;
+ typedef Eigen::VectorXd Vec;
+ typedef Eigen::Map<Vec> MVec;
+ typedef Eigen::ArrayXXd Ar2;
+ typedef Eigen::Map<Ar2> MAr2;
+ typedef Eigen::MatrixXd Mat;
+ typedef Eigen::Map<Mat> MMat;
+ protected:
+ // lme4Eigen::merPredD *d_pred;
+ // lme4Eigen::lmResp *d_resp;
+ MVec d_dev;
+ MMat d_fixef;
+ MVec d_sigma;
+ MMat d_ranef;
+ public:
+ // all the work is done in the constructor
+ mcmcsamp(lme4Eigen::merPredD *pred, lme4Eigen::lmResp *resp,
+ SEXP dev, SEXP fixef, SEXP sigma, SEXP ranef);
+ };
+}
+
+#endif /* LME4_GLMFAMILY_H */
+
Modified: pkg/lme4Eigen/src/predModule.cpp
===================================================================
--- pkg/lme4Eigen/src/predModule.cpp 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/src/predModule.cpp 2012-02-28 23:26:26 UTC (rev 1629)
@@ -13,6 +13,8 @@
using std::invalid_argument;
using std::runtime_error;
+ using Eigen::ArrayXd;
+
typedef Eigen::Map<MatrixXd> MMat;
typedef Eigen::Map<VectorXd> MVec;
typedef Eigen::Map<VectorXi> MiVec;
@@ -263,6 +265,23 @@
std::copy(newU0.data(), newU0.data() + d_q, d_u0.data());
}
+ template <typename T>
+ struct Norm_Rand : std::unary_function<T, T> {
+ const T operator()(const T& x) const {return ::norm_rand();}
+ };
+
+ inline static VectorXd Random_Normal(int size, double sigma) {
+ return ArrayXd(size).unaryExpr(Norm_Rand<double>()) * sigma;
+ }
+
+ void merPredD::MCMC_beta_u(const Scalar& sigma) {
+ VectorXd del2(d_RX.matrixU().solve(Random_Normal(d_p, sigma)));
+ d_delb += del2;
+ VectorXd del1(Random_Normal(d_q, sigma) - d_RZX * del2);
+ d_L.solveInPlace(del1, CHOLMOD_Lt);
+ d_delu += del1;
+ }
+
VectorXi merPredD::Pvec() const {
int* ppt((int*)d_L.factor()->Perm);
VectorXi ans(d_q);
Modified: pkg/lme4Eigen/src/predModule.h
===================================================================
--- pkg/lme4Eigen/src/predModule.h 2012-02-28 14:41:34 UTC (rev 1628)
+++ pkg/lme4Eigen/src/predModule.h 2012-02-28 23:26:26 UTC (rev 1629)
@@ -170,6 +170,7 @@
int info() const {return d_L.info();}
void installPars(const Scalar& f);
+ void MCMC_beta_u(const Scalar& sigma);
void setBeta0(const VectorXd&);
void setTheta(const VectorXd&);
void setU0(const VectorXd&);
More information about the Lme4-commits
mailing list