Rcpp Armadilllo 中的模板 类 arma::Col

Template classes arma::Col in Rcpp Armadilllo

亲爱的程序员们,

为了为大学研讨会编写 bootstrapped 回归版本,我尝试将我的 R 版本实现到 Rcpp 中(这是为了将 R 与 C++/Rcpp 进行比较)。但是,我一直坚持收到的错误消息,因为我并不真正理解这些(C++ 新手,尤其是 Armadillo 的挣扎)。

这是我正在使用的代码。第一个函数是我从互联网上复制的函数,用于获取矩阵的适当行子集(以实现适当的非参数 bootstrap):

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

arma::mat testrows(arma::mat x, arma::Col idx) {
  arma::mat xsub;
  xsub = x.rows(idx-1);
  return xsub;
}

此外,这是我为实际 bootstrap 版本编写的代码:

   #include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat betaBoot(arma::colvec& y, const arma::mat& X, const int nboot) {
    int n = X.n_rows, k = X.n_cols;
    arma::mat betaHat(k,nboot);
    for(int i = 0; i < nboot; i++){
        Rcpp::NumericVector colId = Rcpp::runif(n, 0, (n-1));
        arma::mat X_boot = testrows(X, colId);
        arma::colvec y_boot = y(colId);
        betaHat.col(i) = arma::solve(X_boot, y_boot);
    }
    return betaHat;
}

betaHat 是包含 bootstrapped 系数向量(在每一列中)的矩阵,该矩阵的维度应为 k x nboot。 X_boot(在循环内)应该是 bootstrapped 数据和 y_boot 相应的相关观察。 colId 应该是 bootstrapping 过程的随机索引。最后,应该返回betaHat

附件是我在使用 sourceCpp 时收到的错误图片。

也许这是我看不到的简单的东西,或者可能是缺乏经验,但是学习这些东西会很棒。所以,如果你能帮我解决这个问题,那就太好了。 先感谢您 艾琳

编辑:我的 bootstrapped 回归函数现在看起来(和工作)的方式:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat betaBoot(arma::colvec& y, const arma::mat& X, const int nboot) {
    int n = X.n_rows, k = X.n_cols;
    arma::mat betaHat(k,nboot);
    for(int i = 0; i < nboot; i++){
        arma::uvec colId = arma::randi<arma::uvec>(n, arma::distr_param(0,n-1));
      
        arma::mat X_boot = X.rows(colId);
        arma::colvec y_boot = y(colId);
        betaHat.col(i) = arma::solve(X_boot, y_boot);
    }
    return betaHat;
}

感谢我收到的所有有用评论。

来自评论;函数 testrows 为索引采用 arma::Col,但在 betaBoot 中传递了 Rcpp::NumericVector,因此出错。您还使用实数(随机统一)而不是整数索引进行子集化。 RcppRcppArmadillo 提供了 sample 函数,可以在这里使用。 (同样来自犰狳帮助页面应该是uvec)。

您可以使用 Rcpp 糖函数生成 select 行的索引:

Rcpp::IntegerVector x = Rcpp::seq(0, n-1);
arma::uvec idx = Rcpp::as<arma::uvec>(Rcpp::sample(x, n, true)) ;

或者您可以直接使用 arma::randi:

arma::uvec idx = arma::randi<arma::uvec>(n,  arma::distr_param(0,n-1));

使用第二种方法的一些代码,省略了 testrows 作为“我发现我不需要我的(偷来的)助手 函数测试行":

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
        
// [[Rcpp::export]]
arma::mat betaBoot( arma::mat X, arma::vec y, int nboot){
  
  int n = X.n_rows,  k = X.n_cols;
  arma::mat betaHat(k, nboot);
  
  for(int i=0; i < nboot; i++){
    arma::uvec idx = arma::randi<arma::uvec>(n,  arma::distr_param(0,n-1));
    arma::mat X_boot = X.rows(idx);
    arma::vec y_boot = y.elem(idx);
    betaHat.col(i) = arma::solve(X_boot, y_boot);
  }  
  return betaHat;
}

/***R
set.seed(1)
n = 25
nc = 2
x = cbind(1, matrix(rnorm(n*nc), nc=nc))
y = rnorm(n)

sim = betaBoot(x, y, 2000)
*/