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:

First let me detail the sampler I was working on, it’s a reversible jump MCMC with two dimension changing vectors (which is coded as a NumericVector) and (which is coded as a arma::vec). always contains the entries 0 and the maximum observed event time, so this NumericVector has size greater than 1 always. But 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 and . 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 is**

*if(JA>0) *

**If has dimension 1, we cannot delete anymore!**

*{*

* if(JA==1){ * ** has dim 2, so we need to use our special function**

* lam1=Death1( Y, I, s1, lam1, mu, sig, phi, clam1); * **MCMC for **

* if(lam1.size()==1){ ***If the MCMC deleted an entry, delete and entry of **

* s1.erase(1); ***This erases the middle entry of **

* }*

* } *

* else{*

* Death3( Y, I, s1, lam1, mu, sig, phi, clam1); *

**sampler if**

*}*

*}*

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

Additionally, we have the same issue for adding an entry to the 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 **

* 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 **

……………………….

… 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 vector and our proposed additional split point to .**

*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 and the second two entries return the new 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 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 entries**

*save1[1]=save[2];*

*Extracts 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**

*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 or if the length of the returned vector from the sampler is 2. Otherwise it will extract the elements of the returned NumericVector to update amd .

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

Happy Coding!

-Andrew G. Chapple