Skip to content

Has anyone tried Whit Armstrong's CppMCMC library?

4 messages · Douglas Bates, Whit Armstrong

Original
Douglas Bates
# 
Whit Armstrong has written a C++ library called CppMCMC that provides
BUGS-like functionality in compiled code.  I see that one of his
examples uses the contagious bovine pleuropneumonia data from the cbpp
data set in the lme4 package, fitting an overdispersed binomial-like
model to it.

Does anyone have experience with these classes?  This is not directly
an R question as, at present, I don't know of a bridge between CppMCMC
and R (although it would be an interesting project to use Rcpp for
this).
Whit Armstrong
# 
Thanks for the the post, Professor Bates.

CppBugs can be found here: http://github.com/armstrtw/CppBugs

It's a combination of WinBUGS and PyMC.

The api is still very unstable as I'm attempting to navigate between
good design and maintaining familiar WinBUGS conventions.

Rcpp is the best way to use it from R.  I'll add an example of calling
it this way on the front page so that everyone can try it.

My thinking is that since it's already very inconvenient to run a BUGS
model from R (writing the BUGS script, providing inits, passing data
back and forth), it's not really asking too much for users to write a
small c++ class inline in an R script that is the CppBugs equivalent
of a BUGS model.

I'm very keen to have help on this project, so if anyone is interested
please contact me.

I'll send a follow up post highlighting a few of the features of
CppBugs and examples of WinBUGS vs PyMC conventions that I've followed
in the design.

-Whit
On Tue, Nov 2, 2010 at 11:08 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
1 day later
Whit Armstrong
# 
Below is a minimal example of a simple linear model using CppBugs,
Rcpp, and inline.

The cppbugs dir needs to be on the include path of the compiler for
this example to work.  Otherwise, you will need to change the #include
statements to use the local path of your cppbugs install.

Other dependencies you need are Armadillo and Boost.

feedback welcome.

-Whit

require(inline)
require(Rcpp)

cppbugs.model <- '
#include <armadillo>
#include <cppbugs/cppbugs.hpp>

using namespace arma;
using namespace cppbugs;

class TestModel: public MCModel {
public:
  const mat& y; // given
  const mat& X; // given

  Stochastic<vec> b;
  Stochastic<double> tau_y;
  Deterministic<mat> y_hat;
  Stochastic<mat> likelihood;
  Deterministic<double> rsq;

  TestModel(const mat& y_,const mat& X_): y(y_),
X(X_),b(randn<vec>(X_.n_cols)),tau_y(1),y_hat(X*b.value),likelihood(y_,true),rsq(0)
  {
    add(b);
    add(tau_y);
    add(y_hat);
    add(likelihood);
    add(rsq);
  }

  void update() {
    y_hat.value = X*b.value;
    rsq.value = as_scalar(1 - var(y - y_hat.value) / var(y));
    b.dnorm(0.0, 0.0001);
    tau_y.dunif(0,100);
    likelihood.dnorm(y_hat.value,tau_y.value);
  }
};
'

src <- '
mat X(REAL(XR),::Rf_nrows(XR),::Rf_ncols(XR));
mat y(REAL(yr),::Rf_nrows(yr),::Rf_ncols(yr));
int iterations_ = as<int>(iterations);
int burn_ = as<int>(burn);
int thin_ = as<int>(thin);

TestModel m(y,X);
m.sample(iterations_, burn_, thin_);
return Rcpp::List::create(Rcpp::Named("b", m.b.mean()),
Rcpp::Named("tau.y", m.tau_y.mean()), Rcpp::Named("ar",
m.acceptance_ratio()), Rcpp::Named("rsq", m.rsq.mean()));
'

NR <- 100
NC <- 2
icept <- 10
y = rnorm(NR,icept)
X = matrix(rnorm(NR*NC),ncol=NC)
X[,1] <- 1.0
X[,2] <- y + - icept + rnorm(NR);

bayes.reg <- cxxfunction(signature(XR="numeric",
yr="numeric",iterations="integer",burn="integer",thin="integer"),
body=src, include=cppbugs.model, plugin="Rcpp")
ans <- bayes.reg(X,y,1e5L,1e4L,10L)
print(ans)
On Tue, Nov 2, 2010 at 1:02 PM, Whit Armstrong <armstrong.whit at gmail.com> wrote:
Douglas Bates
# 
One of the first things that I would change is to use the
RcppArmadillo plugin instead of the Rcpp plugin.  That way you can
have the inline package find the armadillo headers for you.  Long term
I would suggest creating a derived package RcppMCMC (or whatever name
seemed appropriate) and providing a plugin to the inline package so
that the headers and libraries do not need to be declared explicitly.

In using the C macros and functions REAL, Rf_nrows, ... you are
leaving yourself open to tromping on memory that doesn't belong to
you.  It is not too bad a risk when using the inline package because
you specify a signature.  However, if you build a more general package
and don't use inline then there is a risk that you will get something
passed that is not a numeric matrix.    I would go through the Rcpp
classes instead because they do all the checking of appropriate types
for you.

That is, I would rewrite the first line as

Rcpp::NumericMatrix x(XR);
mat X(x.begin(), x.nrows(), x.ncols(), false);

The false at the end of the constructor indicates that there is no
need to copy the storage.

And do you really need a matrix for y or is it suitable to use a
column vector?  Consider the following code from one of the examples
in the RcppArmadillo package

	Rcpp::NumericVector yr(ys);			// creates Rcpp vector from SEXP
	Rcpp::NumericMatrix Xr(Xs);			// creates Rcpp matrix from SEXP
	int n = Xr.nrow(), k = Xr.ncol();

	arma::mat X(Xr.begin(), n, k, false);   	// reuses memory and avoids extra copy
	arma::colvec y(yr.begin(), yr.size(), false);

	arma::colvec coef = arma::solve(X, y);      	// fit model y ~ X
	arma::colvec res = y - X*coef;			// residuals
On Wed, Nov 3, 2010 at 1:55 PM, Whit Armstrong <armstrong.whit at gmail.com> wrote: