我遇到了我无法缠住头的事情。这是较大的编码工作的一部分,但此处有一个最小的示例:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List foo(arma::vec & tau2, const arma::vec & nu) {

  arma::vec bet = Rcpp::rnorm(3);
  tau2 = R::rgamma(1, arma::as_scalar(sum(pow(bet, 2)/nu)));

  return Rcpp::List::create(Rcpp::Named("nu") = nu,
                            Rcpp::Named("tau2") = tau2);
}

(tau2虽然是一个标量,但在这里是 vector ,因为我想通过引用传递它:function pass by reference in RcppArmadillo)

让我感到困惑的是,如果我现在运行以下R代码:
n <- 3
m <- matrix(0, n, 1)
for (r in 1:1000) {
  tau2 <- 1.0
  nu <- matrix(1, n, 1)
  upd <- foo(tau2, nu)
}

我得到:
error: element-wise division: incompatible matrix dimensions: 3x1 and 18x1
Error in foo(tau2, nu) :
  element-wise division: incompatible matrix dimensions: 3x1 and 18x1
18x1在哪里变化;大多数情况下是0x1,但始终是3的倍数。

查看输出:
> nu
         [,1]     [,2]     [,3]     [,4]
[1,] 4.165242 4.165242 4.165242 4.165242
[2,] 4.165242 4.165242 4.165242 4.165242
[3,] 4.165242 4.165242 4.165242 4.165242
> upd
$nu
     [,1]
[1,]    1
[2,]    1
[3,]    1

$tau2
         [,1]
[1,] 4.165242

也就是说,尽管将nu声明为常量引用(之所以这样做是因为我不希望更改它),但它还是被更改了。它填充的值是upd$tau2(但是为什么呢?)。

奇怪的是,我可以通过看似无意义的更改来使行为消失:
  • tau2 <- 1.0nu <- matrix(1, n, 1)(或两者)放在循环
  • 之外
  • 删除第一个参数中的引用(即使用arma::vec tau2)
  • 不将pow(bet, 2)除以nu
  • 更改为nu <- rep(1, n)

  • 也许最令人困惑的部分是,如果我在循环内部选择代码块并重复运行它,那么它将起作用(!)。但是,如果我使用循环运行R代码,它将在第二次迭代时崩溃。

    因为我似乎能够解决问题,所以我最有兴趣了解此处的情况。我怀疑这仅仅是由于我缺乏C++专业知识以及各种变量类型的鲁ck性造成的,因此了解导致这些情况的原因将非常有值(value)。

    最佳答案

    两种修复:

  • tau2是 double (这里是@dirkeddelbuettel的模仿)
  • 在保存到NumericVector之前生成的长度为n的bet的临时变量

  • 码:
    #include <RcppArmadillo.h>
    // [[Rcpp::depends(RcppArmadillo)]]
    
    // [[Rcpp::export]]
    Rcpp::List foo(double tau2, const arma::vec & nu) {
    
      int n = nu.n_elem;
    
      Rcpp::NumericVector x = Rcpp::rnorm(n);
      arma::vec bet = arma::vec(x.begin(), n, true, false);
    
      tau2 = R::rgamma(1, arma::as_scalar(sum(pow(bet, 2) / nu)));
    
      return Rcpp::List::create(Rcpp::Named("nu") = nu,
                                Rcpp::Named("tau2") = tau2);
    }
    

    测试用例:
    n <- 3
    m <- matrix(0, n, 1)
    for (r in 1:1000) {
      tau2 <- 1.0
      nu <- matrix(1, n, 1)
      upd <- foo(tau2, nu)
    }
    upd
    #> $nu
    #>      [,1]
    #> [1,]    1
    #> [2,]    1
    #> [3,]    1
    #>
    #> $tau2
    #> [1] 3.292889
    

    关于c++ - 了解Rcpp中的奇怪行为,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/49901812/

    10-12 17:39