[R-sig-ME] Has anyone tried Whit Armstrong's CppMCMC library?

[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:
> 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:
>> 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:
>>> 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).
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>
>