进行平均关系的 Rcpp rank 函数

Question

我有一个从这里偷来的自定义排名函数（经过一些修改）:) Rcpp sugar for rank function

问题是它的关系很少，我需要平均关系

这就是我所拥有的

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector sort_rcpp(NumericVector x) {
  std::vector<double> tmp = Rcpp::as< std::vector<double> > (x);
  std::sort(tmp.begin(), tmp.end());
  return wrap(tmp);
}

// [[Rcpp::export]]
IntegerVector rank_(NumericVector x) {
  return match(x, sort_rcpp(x));
}


/*** R
x <- c(1:5, 1:5)

rank(x, ties = 'average')
rank(x, ties = 'min')
rank_(x)
*/

将其保存到文件然后运行它会得到结果

Rcpp::sourceCpp('~/Documents/Rank.cpp')

返回

# x <- c(1:5, 1:5)
#
# # what I need
# rank(x, ties = 'average')
# [1] 1.5 3.5 5.5 7.5 9.5 1.5 3.5 5.5 7.5 9.5
#
# # What I am getting
# rank(x, ties = 'min')
# [1] 1 3 5 7 9 1 3 5 7 9
#
# rank_(x)
# [1] 1 3 5 7 9 1 3 5 7 9

我需要在 c++ 代码中修改哪些内容以匹配基础 R 中的平均排名函数？

Answer 1

这是 René Richter 在shayaa's link 中的代码的改编版本——主要区别在于使用Rcpp::seq 而不是std::iota 以及处理NA 比较的自定义比较器：

#include <Rcpp.h>

class Comparator {
private:
    const Rcpp::NumericVector& ref;

    bool is_na(double x) const 
    {
        return Rcpp::traits::is_na<REALSXP>(x);    
    }

public:
    Comparator(const Rcpp::NumericVector& ref_)
        : ref(ref_)
    {}

    bool operator()(const int ilhs, const int irhs) const
    {
        double lhs = ref[ilhs], rhs = ref[irhs]; 
        if (is_na(lhs)) return false;
        if (is_na(rhs)) return true;
        return lhs < rhs;
    }
};

// [[Rcpp::export]]
Rcpp::NumericVector avg_rank(Rcpp::NumericVector x)
{
    R_xlen_t sz = x.size();
    Rcpp::IntegerVector w = Rcpp::seq(0, sz - 1);
    std::sort(w.begin(), w.end(), Comparator(x));

    Rcpp::NumericVector r = Rcpp::no_init_vector(sz);
    for (R_xlen_t n, i = 0; i < sz; i += n) {
        n = 1;
        while (i + n < sz && x[w[i]] == x[w[i + n]]) ++n;
        for (R_xlen_t k = 0; k < n; k++) {
            r[w[i + k]] = i + (n + 1) / 2.;
        }
    }

    return r;
}

---

针对base::rank验证结果，

x <- c(1:7, 1:2, 1:5, 1:10)
all.equal(
    rank(x, ties.method = "average"), 
    avg_rank(x)
)
# [1] TRUE 

另请注意，这会正确处理NAs，而您的版本不能：

all.equal(
    rank(c(NA, x), ties.method = "average"), 
    avg_rank(c(NA, x))
)
# [1] TRUE

all.equal(
    rank(c(NA, x), ties.method = "average"), 
    rank_(c(NA, x))
)
# Error: can't subset using a logical vector with NAs

---

这里是上面的向量x 的基准：

microbenchmark::microbenchmark(
    ".Internal" = .Internal(rank(x, length(x), ties = "average")),
    avg_rank(x),
    "base::rank" = rank(x, ties.method = "average"),
    rank_(x),
    times = 1000L
)
# Unit: microseconds
#         expr    min     lq      mean median     uq     max neval
#    .Internal  1.283  1.711  2.029777  1.712  2.139  65.002  1000
#  avg_rank(x)  2.566  3.422  4.057262  3.849  4.277  23.521  1000
#   base::rank 13.685 16.251 18.145440 17.534 18.390 163.360  1000
#     rank_(x) 25.659 28.653 31.203092 29.936 32.074 112.898  1000

这是一个 1e6 长度向量的基准（我没有包括 rank_，因为即使是一次评估也需要永远；见下文）：

set.seed(123)
xx <- sample(x, 1e6, TRUE)

microbenchmark::microbenchmark(
    ".Internal" = .Internal(rank(xx, length(xx), ties = "average")),
    avg_rank(xx),
    "base::rank" = rank(xx, ties.method = "average"),
    times = 100L
)
# Unit: milliseconds
#          expr      min       lq     mean   median       uq      max neval
#     .Internal 302.2488 309.7977 334.7977 322.0396 347.4779 504.1041   100
#  avg_rank(xx) 304.4435 309.9840 330.4902 316.7807 331.6825 427.7171   100
#    base::rank 312.1196 327.3319 351.6237 343.1317 366.7316 445.9004   100

对于更大尺寸的向量，这三个函数的运行时间更接近。理论上，.Internal 调用应该总是比base::rank 快一点，因为它放弃了在后者的主体中进行的额外检查。但是，第二个基准测试中的差异不太明显，因为执行这些检查所需的时间量在函数总运行时间中所占的比例要小得多。

---

附带说明一下，您的代码效率如此之低的一个明显原因是您在循环体中创建向量：

for (int i = 0; i < n; ++i) {
    NumericVector xVal = x[x == x[i]];              // here
    IntegerVector Match = match(xVal, sortX);       // here
    double minM = min(Match);
    int matchSize = Match.size();
    NumericVector aves = NumericVector(matchSize);  // here

    for (int k = 0; k < matchSize; ++k) {
        aves[k] = minM + (k);
    }
    ranks[i] = sum(aves)/Match.size();
}

avg_rank 和 R 的实现（我相信你可以仔细检查the source code）都只创建了两个额外的向量，而不管输入的大小。您的函数正在创建 2 + 3 * N 个向量 (!!!)，其中 N 是您输入的大小。

最后，与性能无关，而是像这样排序（不会正确处理NAs），

NumericVector sort_rcpp(NumericVector x) {
    std::vector<double> tmp = Rcpp::as< std::vector<double> > (x);
    std::sort(tmp.begin(), tmp.end());
    return wrap(tmp);
} 

只需使用 Rcpp 提供的工具：

NumericVector RcppSort(NumericVector x) {
    return Rcpp::clone(x).sort();
}

Answer 2

好的，我在 R 中编写了代码，然后将其翻译为 Rcpp。我希望它和我在问题中的 rank_ 函数（现在是 minrank）一样快，但与基本 R 中的版本相比它相当慢

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
NumericVector sort_rcpp(NumericVector x) {
  std::vector<double> tmp = Rcpp::as< std::vector<double> > (x);
  std::sort(tmp.begin(), tmp.end());
  return wrap(tmp);
}

// [[Rcpp::export]]
IntegerVector minrank(NumericVector x) {
  return match(x, sort_rcpp(x));
}

// [[Rcpp::export]]
NumericVector rank_(NumericVector x) {

  NumericVector sortX = sort_rcpp(x);

  int n = x.size();

  NumericVector ranks = NumericVector(n);

  for(int i = 0; i < n; ++i) {
    NumericVector xVal = x[x == x[i]];
    IntegerVector Match = match(xVal, sortX);
    double minM = min(Match);
    int matchSize = Match.size();
    NumericVector aves = NumericVector(matchSize);

    for(int k = 0; k < matchSize; ++k){
      aves[k] = minM + (k);
    }
    ranks[i] = sum(aves)/Match.size();
  }

  return ranks;

}


/*** R

x <- c(1:7, 1:2, 1:5, 1:10)
r1 <- rank(x, ties = 'average')
r2 <- rank_(x)
all(r1 == r2)

library(microbenchmark)
microbenchmark(
  rank(x, ties = 'average')
  ,rank_(x)
  ,rank(x, ties = 'min')
  ,minrank(x)
  ,.Internal(rank(x, length(x), ties = 'average'))
)

*/


#> x <- c(1:7, 1:2, 1:5, 1:10)
#
#> r1 <- rank(x, ties = 'average')
#
#> r2 <- rank_(x)
#
#> all(r1 == r2)
#[1] TRUE
#
#> library(microbenchmark)
#
#> microbenchmark(
#+   rank(x, ties = 'average')
#+   ,rank_(x)
#+   ,rank(x, ties = 'min')
#+   ,minrank(x)
#+   ,.Internal(rank(x, length(x), ties = 'ave .... [TRUNCATED] 
#Unit: microseconds
#                                            expr    min      lq     mean  median     uq    max neval
#                       rank(x, ties = "average") 13.233 14.6510 17.69987 15.3795 16.432 82.706   100
#                                        rank_(x) 23.035 25.4525 26.98596 26.3180 27.520 42.938   100
#                           rank(x, ties = "min") 13.244 14.3300 17.30062 15.1200 16.763 60.819   100
#                                      minrank(x)  2.529  3.0415  3.66911  3.4265  3.706 14.190   100
# .Internal(rank(x, length(x), ties = "average"))  1.236  1.4185  1.59032  1.4855  1.599  3.125   100
#

我想知道为什么 base::rank 比调用 .Internal rank 慢得多。 Rcpp 中的 minrank 比 base::rank 快得多。

我的 cpp 代码可能非常低效，但它可以工作：|