Rcpp 中的异常行为

Question

请注意，这个错误是从更大的上下文中提取的，显然我不能在这里完全报告。

我在文件中有以下函数fun.cpp

#include <RcppArmadilloExtensions/sample.h>

using namespace Rcpp;

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

arma::vec colMeans(arma::mat data){

   int n_0 = data.n_rows;

   arma::vec xbar(data.n_cols);
   for(int i = 0; i < data.n_rows; i++){
      for(int j = 0; j < data.n_cols; j++){
         xbar[j] += data(i,j) /n_0;
      }
   }

   return xbar;
}

// [[Rcpp::export]]
List PosteriorNIW(arma::mat data, arma::vec mu0, double lambda0,
              double df0, arma::mat V){

   // Compute posterior
   int n = data.n_rows;

   arma::vec xbar = colMeans(data);

   double lambdan = lambda0 + n;

   arma::vec mun = (lambda0 * mu0 + n * xbar) / lambdan;

   arma::mat S;
   S.zeros(data.n_cols, data.n_cols);
   for(int i = 0; i < n; i++){
      S += (arma::conv_to<arma::vec>::from(data.row(i)) - xbar) * arma::trans(arma::conv_to<arma::vec>::from(data.row(i)) - xbar);
   }

   arma::mat Vn = V + S + ((lambda0*n)/(lambda0 + n)) * (xbar - mu0) * arma::trans(xbar - mu0);

   return List::create(_["mun"] = mun,
                  _["Vn"] = Vn,
                  _["lambdan"] = lambdan);
 }

现在打电话：

library(Rcpp); library(RcppArmadillo)
mu0 <- c(3,3)
V0 <- matrix(c(2.5,0.0,0.0,2.5), nrow = 2)
sourceCpp("fun.cpp")

data <- cbind(rep(5,15),rep(0,15))
PosteriorNIW(data,  mu0, 1, 1, V0)

给出了预期的结果。

$mun
    [,1]
[1,] 4.8750
[2,] 0.1875

$Vn
     [,1]    [,2]
[1,]  6.250 -5.6250
[2,] -5.625 10.9375

$lambdan
[1] 16

现在，如果我向文件 fun.cpp 添加以下函数（同样，这些是从更大的上下文中获取的，所以不要费心去理解，只需粘贴它们）奇怪的事情就会发生：

// [[Rcpp::export]]
NumericMatrix myFun(arma::mat t_dish, arma::cube data){
   int l = 0;
   for(int j = 0; j < data.n_rows; j++){
      l++;
   }
   NumericMatrix Dk(l, 2);
   return Dk;
}

// [[Rcpp::export]]
int myFun2(arma::cube n_cust){

   arma::mat temp = n_cust.subcube(arma::span(0), arma::span(), arma::span());
   int i;
   for(i = 0; i < n_cust.n_cols; i++){
      arma::rowvec temp2 = temp.row(i);
   }

   return i + 1;
}

// [[Rcpp::export]]
arma::vec myFun3(arma::mat k_tables){
   arma::vec temp(k_tables.n_cols * k_tables.n_rows);
   int l = 0;
   if(!R_IsNA(k_tables(0,0))){
      l++;
   }

   arma::vec temp2(l);
   arma::vec tmp3 = sort(temp2);
   return tmp3;
}

double myFun4(arma::vec x, double nu, arma::vec mu, arma::mat Sigma){
   arma::vec product = (arma::trans(x - mu) * arma::inv(Sigma) * (x - mu));
   double num = pow(1 + (1 / nu) * product[0], - ( nu + 2 ) / 2);
   double den = pow(sqrt(M_PI * nu),2) * sqrt(arma::det(Sigma));
   return  num / den;
}

bool myFun5(NumericVector X, double z) { 
   return std::find(X.begin(), X.end(), z)!=X.end(); 
}

重复调用 PosteriorNIW(data, mu0, 1, 1, V0) 每次都会给出不同的结果。请注意，函数中没有随机性，显然这些函数没有影响，因为它们没有在原始函数中调用。

我在另一台机器上试过以确保这不是我的编译器的问题，但错误一直在发生。

我知道删除这些功能（即使只是其中一个）可以解决问题，但显然当我使用更多功能时这不是一个可行的解决方案。

我想知道其他用户是否能够复制此行为，如果是，是否有修复程序。

提前致谢

编辑：

R的版本是3.3.2，Rtools是3.4。 Rcpp 和 RcppArmadillo 都是最新的

Answer 1

您没有在 colMeans 函数中对 xbar 进行归零。如果我这样做：

arma::vec colMeans(arma::mat data){

   int n_0 = data.n_rows;

   arma::vec xbar;
   xbar.zeros(data.n_cols);
   for(int i = 0; i < data.n_rows; i++){
      for(int j = 0; j < data.n_cols; j++){
         xbar[j] += data(i,j) /n_0;
      }
   }

   return xbar;
}

我每次都收到这个：

> PosteriorNIW(data,  mu0, 1, 1.1, V0)
$mun
       [,1]
[1,] 4.8750
[2,] 0.1875

$Vn
       [,1]    [,2]
[1,]  6.250 -5.6250
[2,] -5.625 10.9375

$lambdan
[1] 16

即使我确实添加了您的额外代码块。

我不知道这些向量是否被记录为由它们的构造函数初始化为零（在这种情况下这可能是那里的错误），在这种情况下它是你的错误！

Rcpp 中的异常行为

Unexpected behaviour in Rcpp

r

rcpp