执行平均关系的 Rcpp 等级函数
Rcpp rank function that does average ties
我有一个从这里偷来的自定义排名函数（经过一些修改）:)
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')
哪个returns
# 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 的平均排名函数？
好的，我在 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 代码可能效率低得可怕，但它确实有效：|
这是 René Richter 在 中的代码的改编版本 -- 主要区别在于使用 Rcpp::seq 而不是 std::iota 和处理 [=20 的自定义比较器=] 比较：
#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 的实现（我相信，您可以仔细检查 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();
}
我有一个从这里偷来的自定义排名函数（经过一些修改）:) 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')
哪个returns
# 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 的平均排名函数？
好的，我在 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 代码可能效率低得可怕，但它确实有效：|
这是 René Richter 在 Rcpp::seq 而不是 std::iota 和处理 [=20 的自定义比较器=] 比较：
#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 的实现（我相信，您可以仔细检查 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();
}