c++ - 在 Rcpp 中生成整数样本

Question

我想创建一个由 5 个整数组成的随机 vector ，范围为 1:10。我只能使用基本的 Rcpp。 (没有 C 库)

目前我有:

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector test(){
NumericVector z(5);
for (int i=0; i<5 ++i)
z[i] = R::runif(1,10);
return z; 
}
/***R
test()
*/

但是:

• 不是整数

• 它不是唯一的。

Answer 1

这可以用 std::random_shuffle 简洁地完成:

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::IntegerVector sample_int() {
    Rcpp::IntegerVector pool = Rcpp::seq(1, 10);
    std::random_shuffle(pool.begin(), pool.end());
    return pool[Rcpp::Range(0, 4)];
} 

示例输出:

sample_int()
# [1] 9 2 5 1 7

sample_int()
# [1]  1 10  5  3  8

sample_int()
# [1] 5 9 3 2 8 

---

郑重声明，您的代码没有返回整数，因为

• ::runif 返回 double 值；和
• 您的函数的返回类型是 NumericVector 而不是 IntegerVector

---

虽然在处理小范围(例如您的示例中使用的范围 (1, ..., 10))时它是无关紧要的，但这种方法不是很有效(特别是当被采样的元素数量远小于绘图池)，因为 std::random_shuffle 随机播放整个范围。通过几个辅助函数，我们可以做得更好(假设 std::rand 对您的目的而言“足够”随机):

#include <Rcpp.h>

// C++ 98
template <typename Iter, typename T>
inline void iota(Iter first, Iter last, T value) {
    while (first != last) {
        *first++ = value++;
    }
}

template <typename T>
inline T pop_random(std::vector<T>& v) {
    typename std::vector<T>::size_type pos = std::rand() % v.size();
    T res = v[pos];

    std::swap(v[pos], v.back());
    v.pop_back();

    return res;
}

// [[Rcpp::export]]
Rcpp::IntegerVector sample_int2(int n, int min, int max) {
    Rcpp::IntegerVector res(n);
    std::vector<int> pool(max + 1 - min);
    iota(pool.begin(), pool.end(), min);

    for (R_xlen_t i = 0; i < n; i++) {
        res[i] = pop_random(pool);
    }

    return res;
}

并将原解进行泛化比较:

// [[Rcpp::export]]
Rcpp::IntegerVector sample_int(int n, int min, int max) {
    Rcpp::IntegerVector pool = Rcpp::seq(min, max);
    std::random_shuffle(pool.begin(), pool.end());
    return pool[Rcpp::Range(0, n - 1)];
}

---

microbenchmark::microbenchmark(
    "sample_int" = sample_int(100, 1, 1e6),
    "sample_int2" = sample_int2(100, 1, 1e6),
    times = 300L
)
# Unit: milliseconds
#         expr       min        lq      mean    median        uq       max neval
#   sample_int 20.639801 22.417594 23.603727 22.922765 23.735258 35.531140   300
#  sample_int2  1.504872  1.689987  1.789866  1.755937  1.830249  2.863399   300

microbenchmark::microbenchmark(
    "sample_int" = sample_int(1e5, 1, 1e6),
    "sample_int2" = sample_int2(1e5, 1, 1e6),
    times = 300L
)
# Unit: milliseconds
#         expr      min        lq      mean    median        uq       max neval
#   sample_int 21.08035 22.384714 23.295403 22.811011 23.282353 34.068462   300
#  sample_int2  3.37047  3.761608  3.992875  3.945773  4.086605  9.134516   300