Check BIC via RcppArmadillo(Rcpp)

Question

I want check the BIC usage in statistics. My little example, which is saved as check_bic.cpp, is presented as follows:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
List check_bic(const int N = 10, const int p = 20, const double seed=0){
    arma_rng::set_seed(seed); // for reproducibility
    arma::mat Beta = randu(p,N); //randu/randn:random values(uniform and normal distributions)
    arma::vec Bic = randu(N);
    uvec ii = find(Bic == min(Bic)); // may be just one or several elements
    int id = ii(ii.n_elem); // fetch the last one element
    vec behat = Beta.col(id); // fetch the id column of matrix Beta
    List ret;
    ret["Bic"] = Bic;
    ret["ii"] = ii;
    ret["id"] = id;
    ret["Beta"] = Beta;
    ret["behat"] = behat;
    return ret;
}

Then I compile check_bic.cpp in R by

library(Rcpp)
library(RcppArmadillo);
sourceCpp("check_bic.cpp")

and the compilation can pass successfully. However, when I ran

check_bic(10,20,0)

in R, it shows errors as

error: Mat::operator(): index out of bounds
Error in check_bic(10, 20, 0) : Mat::operator(): index out of bounds

I check the .cpp code line by line, and guess the problems probably happen at

uvec ii = find(Bic == min(Bic)); // may be just one or several elements
int id = ii(ii.n_elem); // fetch the last one element

since if uvec ii only has one element, then ii.n_elem may be NaN or something else in Rcpp (while it's ok in Matlab), while I dont konw how to deal with case. Any help?

I search solutions on *stack* website and according to [A simple RcppArmadillo index issue](https://stackoverflow.com/questions/29640048/a-simple-rcpparmadillo-index-issue), I change `int id = ii(ii.n_elem); ` to `int id = ii.tail(1);`, then still have the errors `cannot convert arma::subview_col to int in initialization int id = ii.tail(1); ` — John Stone, Jan 31 '20 at 08:07
C++ counts from 0. If `ii.n_elem` is non-zero, the last element of `ii` is `ii.n_elem-1`. Or you can use [.back()](http://arma.sourceforge.net/docs.html#compat_container_fns) instead, like so: `id = ii.back()` — mtall, Jan 31 '20 at 12:22
@mtall is right -- using `back()` will fix your problem. Could be an answer mtall, no? — duckmayr, Jan 31 '20 at 12:33
Wow, quite helpful answers and save my days! Both `.back()` and `ii.n_elem-1` work! Tks a lot!@mtall And also thank you!@duckmayr — John Stone, Jan 31 '20 at 12:55

Check BIC via RcppArmadillo(Rcpp)

0 Answers0