If you’re doing MCMC with multidimensional parameters, you’re inevitable going to have to use matrix multiplication. In c++ this is much more easily said than done, but with the package RCPP armadillo, it’s actually possible. Why is it so hard in c++? Because in c++, matrices are condensed to longer vectors, row by row so to fill them in we have to be careful about how we index the vector of matrix quantities. Since I have to do this for my MCMC++ purposes, I’m going to show two simple examples which are returning a matrix inverse and returning the quadratic product of a vector and a matrix. You can find a list of all the arma:: functions here (which are used in conjunction with Rcpp to make RcppArmadillo)
http://arma.sourceforge.net/docs.html
But first, load the package RcppArmadillo into R with library(RcppArmadillo) and change the header in your c++ file from #include <Rcpp.h> to #include <RcppArmadillo.h> or else you’ll get an error message since RcppArmadillo already loading in Rcpp. Now for a quick example, add this function fix to your c++ file:
//[[Rcpp::export]]
arma::mat fix(arma::mat x) {
arma::mat y;
y=arma::inv(x);
return(y) ;
}
We have to tell c++ that we want to return a matrix with arma::mat fix() and that we have a matrix input arma::mat x.
arma::mat y; creates a matrix y that we will use to store the inverse of x.
y=arma::inv(x); actually does the matrix inversion and stores it into the matrix y.
And I think by now we know what return(y) does. (If not check out Banging My Head on the Keyboard #1). Notice that we have to invoke the command arma:: before any armadillo commands, much like we needed to use R:: to get random numbers in the previous post. Now let’s test it out in R.
For good measure, I did solve(X) so you can see that it does the same thing as in R.
As a secondary example, I have to return the matrix quadratic product of the difference in two vectors and a matrix, which is involved in the multivariate normal distribution. We want to calculate:
For vectors a,b and a matrix X. Note: This returns a SCALAR value, not a matrix so we need to specify that the function returns a double value. But let’s make it a little more complicated to show more
functionality of the armadillo package, assume that b is created in the function as a vector of some double constant \mu that we will input. The c++ code that accomplishes this is:
//[[Rcpp::export]]
double fix1(arma::mat X,arma::vec a, double mu) {
int J=a.n_rows;
arma::vec y(J);
y.ones();
y=y*mu;
arma::vec b = a-y;
arma::mat z= b.t() * X * b;
return(as_scalar(z));
}
The first thing to note is we are returning a scalar variable with double and the function fix 1 has a matrix input X, a vector input a and a double input mu.
int J=a.n_rows; sets an integer value J equal to the number of elements in the vector a. (ncol would return 1 since a n vector is just a matrix with dimension nx1).
arma::vec y(J); creates a vector of length J.
y.ones(); Fills the vector y with all ones, this works for vectors and matrices alive.
y=y*mu; Multiplies the vector elements by the double constant mu, this also works the same with matrices.
This three lines above are equivalent to rep(mu,J) in R.
arma::vec b = a-y; does the difference we want to use for multiplication and returns a vector b.
arma::mat z= b.t() * X * b; creates the product that we want, we MUST store this as a matrix for the calculations to work properly.
Finally, as_scalar(z) converts a 1X1 matrix into a double valued number.
NOTE: z(1,1) won’t work here!! But in general for a matrix X, X(i,j) will index the matrix entry i,j. This DOESN’T WORK for 1X1 matrices!! Save yourself the struggle!!
Here’s what the output looks like in R, with a comparison to R’s usual matrix manipulation:
As you can see, we can add to this quantity too, so it’s really returning a scalar!
(Ooohs and Ahhs from the crowd erupt)
Ok so I lied about just two examples, let’s do one more that actually uses our fix1 function to get the log density of a multivariate normal random variable. Assume for now that we have a covariance matrix that is scaled by some constant c. For my MCMC++ (has a nice ring to it right?) purposes I have some vector \lambda that has a multivariate normal prior distribution with mean vector of all entries \mu and covariance cS for constants c and a matrix S. The posterior distribution for \lambda does not have a closed form (for my likelihood at least) so we need to do a metropolis hastings step, which includes evaluating the prior density of our proposed \lambda^* vector. Notice that for the multivariate normal distribution, the logarithm can be simplified to
since only these pieces involve \lambda and are needed for our MCMC. So let’s code up this log(density) function so we can use it in our sampling schemes.
//[[Rcpp::export]]
double dmvn(arma::vec lam, double mu, double c, arma::mat S){
double LogD;
arma::mat InvS=inv(S);
LogD = -(.5/c)*fix1(InvS,lam,mu);
return(LogD);
}
As you can see, we’ve already done most of the work!
double dmvn(arma::vec lam, double mu, double c, arma::mat S){} specifies that we have one vector input that we’ll evaluate the density on, two real number and a matrix input S. We specify that we’ll return a real number.
arma::mat InvS=inv(S); takes the inverse of our covariance matrix, which we’ll need for the density.
LogD = -(.5/c)*fix1(InvS,lam,mu); Just runs our covariance matrix through the function we just made and multiplies it by -(.5/c) to get what we saw above. The output looks like:
How nifty!! My hopes are that this will be enough Rcpp Armadillo, it took me a few hours to get it down but those are the big kickers in MCMC++ing.
Happy Coding! Check out some of my non-stats posts if you’re interested!
-Andrew G Chapple
adventuresinstatistics 