Has anyone tried Whit Armstrong's CppMCMC library?

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)

Has anyone tried Whit Armstrong's CppMCMC library?

Thread (4 messages)