Rcpp排名函数，平均联系

Question

我有一个自定义的等级函数，我从这里偷了（有一些修改）:) 等级函数的Rcpp sugar

问题是，它确实关系密切，我需要平均关系

这是我所拥有的

#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

这是shayaa链接中RenéRichter代码的改编版本-主要区别在于使用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 

另请注意，这NA在您的版本不正确的情况下可以正确处理：

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的实现（我相信，你可以仔细检查源代码）只能创建两个额外的载体，不管输入的大小。您的函数正在创建2 + 3 * N个向量（!!!），其中N是您输入的大小。

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

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();
}