这在Rcpp中正确吗？

Question

这在Rcpp中正确吗？

3

我想要比较每一列，并在计算后返回所有结果。我尝试编写代码，但结果不合理。因为如果矩阵中有5列，则结果数量将是5*4/2=10而不是5。我认为问题出在代码中的m上，但我不知道它是否正确。谢谢。

library(Rcpp)
sourceCpp(code='
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h> 

double KS(arma::colvec x, arma::colvec y) {
  int n = x.n_rows;
  arma::colvec w = join_cols(x, y);
  arma::uvec z = arma::sort_index(w);
  w.fill(-1); w.elem( find(z <= n-1) ).ones();
  return max(abs(cumsum(w)))/n;
}
// [[Rcpp::export]]
Rcpp::NumericVector K_S(arma::mat mt) {
  int n = mt.n_cols;
  int m = 1;
  Rcpp::NumericVector results(n);
  for (int i = 0; i < n-1; i++) {
    for (int j = i+1; j < n; j++){
      arma::colvec x=mt.col(i);
      arma::colvec y=mt.col(j);
      results[m] = KS(x, y);
      m ++;
    }
  }
  return results;
}
')

set.seed(1)
mt <- matrix(rnorm(400*5), ncol=5)
result <- K_S(t(mt))

> result
[1] 0.0000 0.1050 0.0675 0.0475 0.0650

- Bertram

这段代码中有超过 n 个循环。 - Stéphane Laurent

1

"结果不合理。" 请更具体地说明。它为什么不合理？您是否收到了错误提示？如果是这样，那么是什么错误？根据您的测试输入，您期望得到什么输出？[如果您正在使用随机测试数据，请提供一个种子以便重现。] - Limey

抱歉，我更新了我的问题。@Limey - Bertram

你声明了一个大小为n的向量results。然而你自己说这个向量应该是大小为n*(n-1)/2。 - Stéphane Laurent

你的问题似乎与排列有关。你可以解决它，而不需要使用K_S部分。 - Dirk Eddelbuettel

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Dirk Eddelbuettel · Accepted Answer

您有一些小错误。在修复它时，我刚刚填充了一个类似于 n by n 矩阵的中间版本--这使得索引错误变得明显。同时返回一个 arma::rowvec 也有助于可能的越界索引错误（默认情况下会报错），但最后您（在这种情况下!!）实际上可以只使用 std::vector 进行扩展。

代码

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

double KS(arma::colvec x, arma::colvec y) {
  int n = x.n_rows;
  arma::colvec w = join_cols(x, y);
  arma::uvec z = arma::sort_index(w);
  w.fill(-1); w.elem( find(z <= n-1) ).ones();
  return max(abs(cumsum(w)))/n;
}
// [[Rcpp::export]]
std::vector<double> K_S(arma::mat mt) {
  int n = mt.n_cols;
  std::vector<double> res;
  for (int i = 0; i < n; i++) {
    for (int j = i+1; j < n; j++){
      arma::colvec x=mt.col(i);
      arma::colvec y=mt.col(j);
      res.push_back(KS(x, y));
    }
  }
  return res;
}

/*** R
set.seed(1)
mt <- matrix(rnorm(400*5), ncol=5)
result <- K_S(mt)
result
*/

输出

> Rcpp::sourceCpp("~/git/stackoverflow/73916783/answer.cpp")
> set.seed(1)
> mt <- matrix(rnorm(400*5), ncol=5)
> result <- K_S(mt)
> result
 [1] 0.1050 0.0675 0.0475 0.0650 0.0500 0.0775 0.0575 0.0500 0.0475 0.0600
>