Efficiently sample a single row of an array at random
I can only comment on the C++ code, as I've never used R nor Rcpp.

For an ignorant reader such as myself, it would have helped to omit the using namespace Rcpp;, so that I could see which names come from that library (I did quite a bit of external reading to even understand the code).

Assuming I'm right that sample() returns an array of length k, then we can use range-based for to iterate through its values (since IntegerVector has suitable begin() and end()):

// sample k values from x, without replacement
const Rcpp::IntegerVector ind_index = Rcpp::RcppArmadillo::sample(x, k, false);

// how many different values do they index in pop_ptr?
std::set&lt;int&gt; distinct_container;
for (int i: ind_index) {
    distinct_container.insert(pop_ptr[i]);
}


In accumulate() (which name I dislike, because it sounds like std::accumulate()), we can reduce the scope of i, j, k and pop_ptr quite simply:

for (int i = 0; i &lt; K; ++i) {
    for (int k = 0;  k &lt; num_specs;  ++k) {
        for (int j = 0;  j &lt; perms;  ++j) {
            const int *const pop_ptr = &amp;(pop(sample_one(perms), 0, sample_one(K)));
            res(j, k, i) = sample_n_distinct(specs_C, k + 1, pop_ptr);
        }
    }
}


Also, I'd be inclined to show that the creation of specs_C is more than just a conversion from specs - 1 - that's actually a constructor argument specifying the length of the new vector:

IntegerVector specs_C(specs - 1);


I don't see any code that ever modifies this vector; if it's supposed to remain full of zeros, then it can be declared const.



With those changes, and an automated re-indent, we get something like this for those two functions (I've compiled, but not tested this):

int sample_one(int n)
{
    return n * R::unif_rand();
}

int sample_n_distinct(const Rcpp::IntegerVector&amp; x,
                      int k,
                      const int * pop_ptr)
{
    // sample k values from x, without replacement
    const auto ind_index = Rcpp::RcppArmadillo::sample(x, k, false);

    // how many different values do they index in pop_ptr?
    std::set&lt;int&gt; distinct_container;
    for (int i: ind_index) {
        distinct_container.insert(pop_ptr[i]);
    }

    return distinct_container.size();
}

// [[Rcpp::export]]
arma::Cube&lt;int&gt; accumulate(const arma::Cube&lt;int&gt;&amp; pop,
                           const Rcpp::IntegerVector&amp; specs,
                           int perms,
                           int K)
{
    auto const num_specs = specs.size();
    arma::Cube&lt;int&gt; res(perms, num_specs, K);

    Rcpp::IntegerVector specs_C(specs - 1);

    for (int i = 0;  i &lt; K;  ++i) {
        for (int k = 0;  k &lt; num_specs; ++k) {
            for (int j = 0;  j &lt; perms; ++j) {
                const auto&amp; sampled = pop(sample_one(perms), 0, sample_one(K));
                res(j, k, i) = sample_n_distinct(specs_C, k + 1, &amp;sampled);
            }
        }
    }
    return res;
}




Sadly, I don't have the domain knowledge to suggest more meaningful improvements other than these fairly mechanical changes that improve readability.

In particular, I don't understand why we have to pass specs_C to accumulate() instead of sampling directly from pop_ptr.

I'm also very wary of taking the address of the element returned by pop() - that seems to assume a memory layout that could vary between platforms, for instance.  I'd expect to use one of the subcube view operations to construct the population:

int sample_n_distinct(const Rcpp::IntegerVector&amp; population, int k)
{
    // sample k values from population, without replacement
    const auto sampled = Rcpp::RcppArmadillo::sample(population, k, false);

    // how many different values do have?
    std::set&lt;int&gt; distinct(sampled.begin(), sampled.end());

    return distinct.size();
}

// [[Rcpp::export]]
arma::Cube&lt;int&gt; accumulate(const arma::Cube&lt;int&gt;&amp; pop,
                           const Rcpp::IntegerVector&amp; specs,
                           int perms,
                           int K)
{
    auto const num_specs = specs.size();
    arma::Cube&lt;int&gt; res(perms, num_specs, K);

    for (int i = 0;  i &lt; K;  ++i) {
        for (int k = 0;  k &lt; num_specs; ++k) {
            for (int j = 0;  j &lt; perms; ++j) {
                const auto&amp; sampled = pop.tube(sample_one(perms), sample_one(K));
                const auto sampled_vector
                    = Rcpp::IntegerVector(sampled.begin(), sampled.end());
                res(j, k, i) = sample_n_distinct(sampled_vector, k + 1);
            }
        }
    }
    return res;
}


I may be wrong on exactly what subcube view is required, so pop.tube() may need to be a different call, but you get the idea.  See how much simpler life is when you can use iterators instead of grubbing around with pointers?