Banging My Head Against the Keyboard: Solving Dimension Changing issues in RCPPArmadillo

If you’ve been reading my Banging My Head Against the Keyboard Series, you know I’ve been struggling to use c++ as a R user (using RCPP). The good news is that I’ve successfully coded up my sampler in c++! It took a lot of effort to get this to work, which I’ll explain along with the solution to the problem.

You can find my other tutorials on RCPP here:

  1. Starting RCPP
  2. Basic c++ for MCMC 
  3. More c++ MCMC
  4. RCPPArmadillo Basics
  5. Debugging in RCPP

First let me detail the sampler I was working on, it’s a reversible jump MCMC with two dimension changing vectors \mathbf{s} (which is coded as a NumericVector) and \mathbf{\lambda} (which is coded as a arma::vec). \mathbf{s} always contains the entries 0 and the maximum observed event time, so this NumericVector has size greater than 1 always. But \mathbf{\lambda} can become a 1 dimensional vector, this is where the issue arises in Rcpp armadillo. Let me explain the function I am using to “delete” one entry of \mathbf{s} and \mathbf{\lambda}. Check out the wrapper below:

void Death3(NumericVector Y, NumericVector I, NumericVector &s1, arma::vec &lam1,
double mu, double sig, double phi, double clam1){

}

So I am passing s1 and lam1 by reference with the & symbol and changing these two vectors within this function. This is called during the MCMC at each iteration. However, when we want to make lam1 go from dimension 2 to dimension 1, RcppArmadillo can’t handle this. To get around this issue, I did the following in my MCMC: first I created a function called Death1 which returns a arma:: vec.

arma::vec Death1(NumericVector Y, NumericVector I, NumericVector &s1, arma::vec &lam1,
double mu, double sig, double phi, double clam1){

arma::vec lam(1);

……………………….

… mcmc details ….

……………………….
U=log(U);
if(U<alpha)
{
arma::vec lam(1);      Sets up a new arma::vec
lam[0]=newLam1;    Sets single value of 1 dimensional arma::vec
return(lam);       returns 1-dimensional arma::vec
}
else{
return(lam1);        returns 2-dimensional previous value
}
}

Notice that even though we are passing s1 and lam1 via reference (to save memory), they are NOT being changed in this function. Either the arma::vec lam or arma::vec lam1 are being returned, depending on whether the proposed move is deleted. But what about s1? We aren’t returning it because if we know it contains (0,MAX), then all we have to do is delete the middle entry of s1=(0,deleted,MAX). We account for this in the MCMC in the following way:

J1=lam1.n_rows-1;      First we want to check how long \mathbf{\lambda} is

if(JA>0)           If \mathbf{\lambda} has dimension 1, we cannot delete anymore!
{

if(JA==1){        \mathbf{\lambda} has dim 2, so we need to use our special function
lam1=Death1( Y, I, s1, lam1, mu, sig, phi, clam1);        MCMC for \mathbf{\lambda}
if(lam1.size()==1){     If the MCMC deleted an entry, delete and entry of \mathbf{s}
s1.erase(1);   This erases the middle entry of \mathbf{s}
}
}
else{
Death3( Y, I, s1, lam1, mu, sig, phi, clam1);    sampler if dim(\mathbf{\lambda}) >2
}
}

So we just need to do a “trick” when we are going to propose deleting an entry of \mathbf{\lambda} when it’s dimension is 2.

Additionally, we have the same issue for adding an entry to the \mathbf{\lambda} vector when it has dimension 1, so we have to do some finagling with c++. The following function accomplishes this:

//[[Rcpp::export]]

NumericVector Birth1(NumericVector Y, NumericVector I, NumericVector &s1, arma::vec &lam1,
double mu, double sig, double phi, double clam1){

double Birth = R::runif(0,s1[1]);      Proposes adding a entry to \mathbf{s}
NumericVector s = s1;
s.push_back(Birth);   Adds the proposal to our vector
std::sort(s.begin(),s.end());   Sorts the vector

arma::vec l2(2);    This sets up our proposed addition to \mathbf{\lambda}

……………………….

… mcmc details ….

……………………….
double U = R::runif(0,1);
U=log(U);
int dec=U<alpha;
if(dec==1)
{
NumericVector save(3);     This creates a NumericVector that is going to contain the \mathbf{\lambda} vector and our proposed additional split point to \mathbf{s}.

save[0]=Birth;   Sets first entry = to the proposed addition
save[1]=l2[0];   The next two stores our new lambda vector.
save[2]=l2[1];

return(save);

}
else{
NumericVector save2(2);
save2[0]=1;     These don’t matter, we just need to return some NumericVector so c++ doesn’t crash.
save2[1]= 1;
return(save2);
}
}

So what this does in a nutshell is return one of two NumericVectors. The first NumericVector (if the proposal is accepted) returns the new entry that we will add to \mathbf{s} and the second two entries return the new \mathbf{\lambda} vector, so we will have to extract these individually. If the move is rejected, we return a junk NumericVector, which has length 2 instead of 3. We will use the length of the vector returned as a proxy for determining whether the parameter vectors should be updated or not.

We use this function in our MCMC in the following manner:

if(J1<Jmax)
{
if(J1==0){   If the dimension of \mathbf{\lambda} is 1, use special function.

save=Birth1( Y, I, s1, lam1, mu, sig, phi, clam1 );     Saves our 2 or 3 dim NumericVector
if(save.size()>2){    If the size is bigger than 2, the birth move was accepted
save1[0]=save[1];     Extracts \mathbf{\lambda} entries
save1[1]=save[2];     Extracts \mathbf{\lambda} entries 
lam1=fillvec(save1);   This is a separate function I made that fills arma::vec with numeric Vectors
s1.push_back(save[0]);      Adds the birth point to \mathbf{s}
std::sort(s1.begin(),s1.end());  Sorts the new vector
}

}
else{
Birth3( Y, I, s1, lam1, mu, sig, phi, clam1 );    Other birth function
}
}

So this function doesn’t change \mathbf{\lambda} or \mathbf{s} if the length of the returned vector from the sampler is 2. Otherwise it will extract the elements of the returned NumericVector to update  \mathbf{\lambda} amd \mathbf{s}.

This sorted through my issues with RcppArmadillo! Let me know if you have difficulty following along!

Happy Coding!

-Andrew G. Chapple

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s