Rcpp::rmultinom 带有指向矩阵子视图的指针

Question

要从 Rcpp 中的多项式分布中抽样，我们可以这样做：

int n = 100;
int k = 3;
arma::vec probs = {0.4,0.2,0.4}
arma::irowvec c(k);
Rcpp::rmultinom(n, probs.begin(), k, c.begin());

我想知道当 C 是矩阵时我们是否可以做同样的事情。我试试

int n = 100;
int k = 3;
arma::vec probs = {0.4,0.2,0.4}
# C (arma::mat C) passed by reference to the function
Rcpp::rmultinom(n, probs.begin(), k, C.row(1).begin());

但是它抛出一个错误。有简单的方法吗？

我想尝试第二种方法，因为我有一个大矩阵 C，我通过引用传递给我的函数，然后我想在多项式之后更新它的行。

MWE：

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(cpp11)]]

using namespace Rcpp;
using namespace arma;

void subfunction(const arma::imat& C) {

  int n = 100;
  int k = 3;
  arma::vec probs = {0.4,0.2,0.4};

  rmultinom(n, probs.begin(), k, C.col(1).begin());
}

// [[Rcpp::export]]
arma::imat myfunction(){

  arma::imat C = { {1, 2}, 
                  {3, 4},
                  {5, 6}};
  subfunction(C);
  Rcpp::Rcout << "C: " <<  C << std::endl;

}

错误是：

test_multim.cpp:18:44: error: ‘class arma::subview_col<int>’ 
has no member named ‘begin’
rmultinom(n, probs.begin(), k, Ct.col(1).begin());

Answer 1

这个问题有两个问题：

rmultinom 函数定义所需的类型。
正在访问 armadillo 矩阵的内存指针。

首先，请注意其中一个错误是：

mondayso.cpp:16:3: error: no matching function for call to 'Rf_rmultinom'
  rmultinom(1, probs.begin(), k, C.colptr(1));
  ^~~~~~~~~
/Library/Frameworks/R.framework/Resources/include/Rmath.h:468:6: note: candidate function not viable: no known conversion from 'double *' to 'int *' for 4th argument
void    rmultinom(int, double*, int, int*);
        ^
1 error generated.

本质上，rmultinorm 函数必须在第 4 个参数上传递一个整数。由于 arma::mat being double by default, the matrix's type is inappropriate. In this case, the C matrix must be arma::imat since it uses the arma::sword 或签名 int 组件的构建。

接下来是armadillo matrices is stored in a column-by-column order (See Wikipedia's entry for details). This means that pointers can only be easily established by column via .colptr的数据。这解决了出现的第二个错误：

error: no member named 'begin' in 'arma::subview_row<int>'
  rmultinom(n, probs.begin(), k, C.row(1).begin());
                                 ~~~~~~~~ ^
1 error generated.

话虽如此，我构建了一个有助于转换的示例。

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
arma::imat test() {

  int n = 100;
  int k = 3;
  arma::vec probs = {0.4,0.2,0.4};

  arma::imat C = { {1, 3, 5},
                  {2, 4, 6} };

  arma::imat Ct = C.t();

  // C++ indices start at 0 (thus, this is the second column!)
  rmultinom(n, probs.begin(), k, Ct.colptr(1));

  return Ct;
}

测试

set.seed(111)
test()
#      [,1] [,2]
# [1,]    1   43
# [2,]    3   18
# [3,]    5   39

Rcpp::rmultinom 带有指向矩阵子视图的指针

Rcpp::rmultinom with a pointer to a matrix subview

armadillo

rcpp