In many analyses, reproducibility is just as important as the quality of the analytical methods used. The ability to reproduce the same set of numbers is a core mechanism of R.
random_stream exists to allow that reproducibility from
R, but entirely in C++.
Creating a stream
The default constructor seeds a new stream from R’s current seed state:
The consequence is the one that matters for users: a function that builds a stream is reproducible from R.
set.seed(1)
one_unif()
#> [1] 0.8642478
set.seed(1)
one_unif() # Reproduces the same number under the same R seed
#> [1] 0.8642478Without setting a seed, we get random numbers as expected.
replicate(10, one_unif())
#> [1] 0.35822373 0.24848048 0.07650453 0.86536336 0.28446347 0.09335021
#> [7] 0.11783802 0.77373274 0.63778202 0.97516822Because constructing a stream consumes R draws, two fresh
streams under one set.seed() are different - exactly as you
would want, and exactly like calling runif() twice:
Reusing one stream
When doing repeated sampling, it’s more efficient to generate the
random_stream once and draw from it many times, rather than
constructing a new one per draw.
[[cppally::register]]
r_vector<r_dbl> unif_deviates(int n){
random_stream r; // <-- Generated once here
r_vector<r_dbl> out(n);
for (r_size_t i = 0; i < n; ++i){
out.set(i, r.unif());
}
return out;
}
set.seed(1)
unif_deviates(10)
#> [1] 0.8642478 0.3523587 0.2399782 0.3767446 0.8174655 0.9434300 0.9112013
#> [8] 0.4336280 0.2818099 0.0432162unif() returns a double in
[0, 1) by default, but can return any range in
[a, b).
Explicit seeds
A stream can also be built from a seed you choose, in which case R’s state is irrelevant and never touched:
[[cppally::register]]
r_vector<r_dbl> unif_seeded(int n, uint64_t seed){
random_stream r(seed);
r_vector<r_dbl> out(n);
for (r_size_t i = 0; i < n; ++i){
out.set(i, r.unif());
}
return out;
}
set.seed(1)
unif_seeded(5, seed = 42)
#> [1] 0.8143051 0.3188210 0.9838942 0.7011356 0.7935045
set.seed(99)
unif_seeded(5, seed = 42) # Same numbers, different R seed
#> [1] 0.8143051 0.3188210 0.9838942 0.7011356 0.7935045Every stream can return the seed it was created with, so a run seeded from R can still be replayed later:
index()
random_stream::index() samples a random integer between
the specified bounds (inclusive at both ends).
set.seed(1)
random_index(0, 1000)
#> [1] 865Unlike the <random> distributions discussed below,
index() is platform-independent. It uses Lemire’s nearly
divisionless method over pure integer arithmetic, so the same seed gives
the same integers everywhere.
Example: efficient sample with replacement
[[cppally::register]]
r_vector<r_int> cpp_sample_int_with_replacement(int n, int size){
random_stream r;
r_vector<r_int> out(size);
for (r_size_t i = 0; i < size; ++i){
auto idx = r.index(0, n - 1);
out.set(i, idx);
}
return out;
}Benchmark comparison against
sample.int(..., replace = TRUE):
library(bench)
mark(
cppally_sample_int_with_replacement = cpp_sample_int_with_replacement(10^6, 10^5),
base_sample_int_with_replacement = sample.int(10^6, 10^5, replace = TRUE),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 cppally_sample_int_with_replac… 224.32µs 370.54µs 2686. 391KB 19.0
#> 2 base_sample_int_with_replaceme… 2.31ms 2.33ms 425. 391KB 4.13In this simple benchmark we achieve a large speed improvement over
base R. The whole loop runs without re-entering R’s RNG once - only the
random_stream r; line does that.
Using <random> distributions
random_stream satisfies
std::uniform_random_bit_generator, so it can be handed
straight to any C++ distribution or algorithm -
std::shuffle, std::normal_distribution and so
on.
#include <random>
[[cppally::register]]
r_vector<r_dbl> normal_deviates(int n){
random_stream r;
r_vector<r_dbl> out(n);
std::normal_distribution<> dist;
for (r_size_t i = 0; i < n; ++i){
out.set(i, dist(r));
}
return out;
}
set.seed(1)
normal_deviates(10)
#> [1] -0.36860578 0.90939231 -0.63676191 -1.34332485 -0.13620834 0.84386574
#> [7] 0.05412269 0.61586757 1.17028316 -1.23532499Two caveats are worth knowing before you rely on this.
The C++ standard specifies the interface of the
<random> distributions but not their algorithms, so
results are not guaranteed to be identical across platforms, compilers,
or standard library versions. The same seed reproduces the same numbers
on one machine, but not necessarily on somebody else’s. Relatedly,
Writing R Extensions asks compiled code not to depend on the
C++11 random number library - which is why random_stream’s
own methods (unif(), index()) are implemented
directly rather than delegating to <random>.
If cross-platform reproducibility matters, prefer
unif()/index(), or R’s own distribution
functions via draw_from_r().
Calling R C API RNG functions
⚠️ To call R’s C API RNG functions - unif_rand(),
norm_rand(), R_unif_index(), or anything from
Rmath.h - wrap them in draw_from_r().
[[cppally::register]]
r_vector<r_dbl> r_normal_deviates(int n){
r_vector<r_dbl> out(n);
draw_from_r([&]{
for (r_size_t i = 0; i < n; ++i){
double v = norm_rand();
out.set(i, v);
}
});
return out;
}
set.seed(1)
r_normal_deviates(5)
#> [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
set.seed(1)
rnorm(5) # Identical - these are R's own draws
#> [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078draw_from_r() loads R’s RNG state before running your
callable and writes it back after, i.e. it is an RAII wrapper around
GetRNGstate() and PutRNGstate(). Without it,
draws never reach .Random.seed, so R hands the same values
out again on its next call.
It returns whatever the callable returns, so it wraps expressions as well as statements:
Wrap the whole loop, as above, rather than each individual draw - the state round trip is the expensive part.
Bootstrapping
Given that we can generate random indices very efficiently, it stands to reason that we can also generate bootstrap samples efficiently.
[[cppally::register]]
r_vector<r_dbl> boot_mean(r_vector<r_dbl> x, int n_boot){
r_size_t n = x.length();
random_stream r;
r_vector<r_dbl> out(n_boot);
for (int b = 0; b < n_boot; ++b){
r_dbl total(0);
// calculate mean
for (r_size_t i = 0; i < n; ++i){
total += x.get(r.index(r_size_t(0), n - 1));
}
out.set(b, total / n);
}
return out;
}boot_mean resamples x (with replacement),
generating n_boot mean estimates. We can then use
quantile to create a percentile confidence interval.
set.seed(42)
x <- rnorm(1000)
set.seed(1)
replicates <- boot_mean(x, 2000)
quantile(replicates, c(0.025, 0.975)) # Percentile confidence interval
#> 2.5% 97.5%
#> -0.08765848 0.03722926Comparison to R - benchmark
We could compare boot_mean to R’s equivalent
mean + sample, but that wouldn’t be a fair
comparison, given that boot_mean is a hand-tuned C++
function which doesn’t generate intermediate vectors.
Instead let’s compare apples to apples and write our own C++ sample-with-replacement.
[[cppally::register]]
r_vector<r_dbl> resample(r_vector<r_dbl> x){
r_size_t n = x.length();
r_vector<r_dbl> out(n);
random_stream r;
for (r_size_t i = 0; i < n; ++i){
out.set(i, x.get(r.index(r_size_t(0), n - 1)));
}
return out;
}Now let’s test it against R.
boot_mean_r <- function(x, n_boot){
vapply(
seq_len(n_boot),
\(i) mean(sample(x, replace = TRUE)),
0
)
}
boot_mean_cppally <- function(x, n_boot){
vapply(
seq_len(n_boot),
\(i) mean(resample(x)),
0
)
}
mark(
cppally = boot_mean_cppally(x, 2000),
base_r = boot_mean_r(x, 2000),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 cppally 23ms 23.1ms 41.8 15.4MB 27.8
#> 2 base_r 63.5ms 63.5ms 15.8 23.1MB 110.The hand-tuned bootstrap mean naturally is faster.
mark(
boot_mean_cppally(x, 2000),
boot_mean(x, 2000),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 boot_mean_cppally(x, 2000) 23.01ms 23.29ms 42.7 15.4MB 28.5
#> 2 boot_mean(x, 2000) 4.36ms 4.36ms 229. 21.8KB 0A flexible and fast bootstrapper
The gap between boot_mean_cppally and
boot_mean is worth pulling apart. Three costs are paid on
every one of those 2,000 iterations: resample() builds a
fresh random_stream, so R’s .Random.seed is
loaded and written back each time - exactly what Reusing one
stream set out to avoid; a new n-length vector is
allocated only to have a single number taken from it; and each replicate
makes a round trip out to R and back.
boot_mean pays none of them per replicate. It constructs
one stream, allocates nothing inside the loop, and never leaves C++.
All three go away together if we hold one stream for the whole run, re-use one intermediate vector, and keep the replicate loop in C++. If we can manage that without fixing the statistic in advance, we get something generic that still performs like the hand-tuned version.
Bootstrap functional:
template <RScalar T, typename FUN>
requires RScalar<std::invoke_result_t<FUN&, const r_vector<T>&>> // Accept only fns that return scalars
auto apply_boot(const r_vector<T>& data, int n_boot, FUN fn){
r_size_t n = data.length();
// Function return type
using stat_t = std::invoke_result_t<FUN&, const r_vector<T>&>;
random_stream r;
r_vector<stat_t> out(n_boot);
// Reused across every replicate
r_vector<T> samples(n);
for (int b = 0; b < n_boot; ++b){
// Sample x using a random index [0, n - 1]
// Set the samples directly into our existing samples vector
for (r_size_t i = 0; i < n; ++i){
auto idx = r.index(r_size_t(0), n - 1);
// data.view() is safe here because we are immediately protecting the result from gc
// by placing it into samples, which is itself protected
samples.set(i, data.view(idx));
}
// Apply function
out.set(b, fn(samples));
}
return out;
}Most of the complexity sits in the signature. The three arguments are:
-
data- the vector being bootstrapped -
n_boot- the number of replicates -
fn- a C++ lambda returning the test-statistic of your choice
The requires clause rejects any fn that
does not return a scalar, and stat_t picks up whatever type
it does return.
Using cppally’s cppally::mean(), we can easily re-write
a new mean bootstrap function.
[[cppally::register]]
r_vector<r_dbl> boot_mean2(r_vector<r_dbl> x, int n_boot){
return apply_boot(x, n_boot, [](const auto& vec){
return mean(vec);
});
}Comparing our new bootstrap mean function against the hand-tuned version.
mark(
boot_mean(x, 2000),
boot_mean2(x, 2000),
check = FALSE
)
#> # A tibble: 2 × 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 boot_mean(x, 2000) 4.35ms 4.36ms 229. 15.7KB 0
#> 2 boot_mean2(x, 2000) 4.45ms 4.46ms 223. 23.5KB 0Almost as fast as the hand-tuned version, which is a nice result
given that boot_mean2 was written using a functional.
One benefit of a generic bootstrap functional like
apply_boot is the ability to quickly write new bootstrap
methods.
// Bootstrapped sum
[[cppally::register]]
r_vector<r_dbl> boot_sum(r_vector<r_dbl> x, int n_boot){
return apply_boot(x, n_boot, [](const auto& vec){
return as<r_dbl>(sum(vec));
});
}
// Bootstrapped count of positive values
[[cppally::register]]
r_vector<r_int> boot_n_positive(r_vector<r_dbl> x, int n_boot){
return apply_boot(x, n_boot, [](const auto& vec){
return vec.reduce([](auto acc, auto curr){
return as<r_int>(acc + (curr > 0));
},
/*init = */ r_int(0));
});
}Note that these two return different types. boot_sum
produces doubles, whereas boot_n_positive counts, and so
produces integers. Neither had to declare that - stat_t
follows whatever the lambda returns, and the output vector is built from
it.
set.seed(42)
boot_sum(x, 5)
#> [1] -71.1128232 18.2339024 0.1102194 -41.3509977 -55.1178594
boot_n_positive(x, 5)
#> [1] 479 487 454 486 479