Rcpp Armadilllo 中的模板类 arma::Col

Posted 2023-02-16

【中文标题】Rcpp Armadilllo 中的模板类 arma::Col

【英文标题】：Template classes arma::Col in Rcpp Armadilllo

【发布时间】：2020-12-24 02:08:21

【问题描述】：

亲爱的程序员们，

为了为大学研讨会编写自举回归版本，我尝试将我的 R 版本实现为 Rcpp（这是为了比较 R 与 C++/Rcpp）。但是，我对收到的错误消息感到困惑，因为我并不真正理解那些（刚接触 C++ 尤其是犰狳的挣扎）。

这是我正在使用的代码。第一个函数是我从互联网上复制来获取矩阵的适当行子集的函数（以实现适当的非参数引导程序）：

    #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;

此外，这是我为实际引导版本编写的代码：

   #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 是一个包含自举系数向量（在每一列中）的矩阵，该矩阵的维度应为 k x nboot。 X_boot（在循环内）应该是自举数据，y_boot 应该是相应的相关观察。 colId 应该是引导过程的随机索引。最后，应该返回betaHat。

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

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

编辑：我的自举回归函数现在看起来（和工作）的方式：

#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;

感谢我收到的所有有用的 cmets。

【问题讨论】：

快速查看testrows 为索引采用arma::Col，但您在betaBoot 中传递了Rcpp::NumericVector。您似乎还使用真实（随机统一）而不是整数索引进行子集（它们应该是 uvec 吗？）

Rcpp 和 RcppArmadillo 都有示例函数

作为一般规则，我建议从最简单的工作构建块开始。在这里，您显然对testrows 有问题，所以我建议您在继续之前能够调用它来随机选择行。

Re: Dirk 建议以testrows() 开头，我认为您应该能够通过not 使用arma::Col 来解决这个问题。如果您查看 excellent Armadillo 文档here，您会发现该规范需要arma::Col<double> 之类的东西；你可能想要arma::vec（相当于arma::Col<double>）

非常感谢大家。事实上，我发现我根本不需要我的（被盗的）辅助函数 testrows，但是我也误解了 .rows 的工作方式（我认为它的工作方式类似于 R 中的 A[n:p,]，而不是 A[c(. ..),])。现在，我直接实现了建议的采样和 .rows 函数。如果你写一个答案，我肯定会接受。再次感谢所有的duckmayr、Dirk 和@user20650。

【参考方案1】：

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

您可以使用Rcpp 糖函数生成索引以选择行：

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