了解 Rcpp 中的奇怪行为

Question

我遇到了一些我无法解决的问题。这是更大的编码工作的一部分，但这里有一个最小的例子：

#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虽然是标量，但这里是向量，因为我想通过引用传递：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.0 或nu <- matrix(1, n, 1)（或两者）置于循环之外
• 删除第一个参数中的引用（即使用arma::vec tau2）
• 不将pow(bet, 2) 除以nu
• 更改为nu <- rep(1, n)

也许最令人困惑的部分是，如果我选择循环内的代码块并重复运行它，它就可以工作（！）。但是，如果我使用循环运行 R 代码，它会在第二次迭代时崩溃。

因为我似乎能够解决问题，所以我最感兴趣的是了解这里发生了什么。我怀疑这只是我缺乏 C++ 专业知识以及对各种变量类型鲁莽的结果，因此了解导致所有这些的原因将非常有价值。

Answer 1

两个修复：

1. tau2 是双精度（此处模仿 @dirkeddelbuettel）
2. 在保存到 bet 之前生成长度为 n 的 NumericVector 的临时变量

代码：

#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

Answer 2

如果我将界面更改为使用double，则一切正常：

代码

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List foo(double & 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);
}

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

演示

R> sourceCpp("/tmp/hejseb.cpp")

R> n <- 3

R> m <- matrix(0, n, 1)

R> for (r in 1:1000) {
+   tau2 <- 1.0
+   nu <- matrix(1, n, 1)
+   upd <- foo(tau2, nu)
+ }
R> upd
$nu
     [,1]
[1,]    1
[2,]    1
[3,]    1

$tau2
[1] 1.77314

R> 

我不确定这些是否是您预期的数字。我真的没有时间完成你想​​要做的事情。