Recreating rmultinom and rpois with Rcpp
This article is originally published at https://gallery.rcpp.org
Sometimes one needs to mimic the exact behavior of R’s Distributions
within C++
code. The incredible Rcpp team has provided access to these distributions through
Rmath.h
(in the R::
namespace), as well as through the Rcpp::
namespace where there
can be two forms: scalar as in R, and vectorized via Rcpp sugar. The behavior of these
functions may not always exactly match what the user expects from the standard R behavior,
particularly if attempting to use the functions in Rmath.h
. In particular, the functions
in Rmath.h
are not vectorized. In what follows, I will use Rcpp to mimic the
behavior of both the rmultinom
and rpois
functions available in base R so that this
functionality and behavior is provided in native C++.
The multinomial distribution
generalizes the binomial distribution to k
discrete outcomes instead of 2; consequently,
it is parameterized in terms of k
probabilities that must sum to 1. The base R function
rmultinom
used for generating multinomial data takes three arguments: n
the number of
simulated data sets to produce, size
, the number of multinomial outcomes to sample for
each data set, and prob
a numeric vector of probabilities. The function returns a k
$\times$ n
integer matrix.
The following C++ code uses the R::rmultinom
function available in Rmath.h
to generate
size
multinomial outcomes. The R::rmultinom
function relies on referencing a pointer
to an IntegerVector
to store the results. We create a helper function, rmultinom_1
,
that draws size
multinomial outcomes from the multinomial distribution based on the
probabilities in prob
. We then do this n
independent times in the function
rmultinom_rcpp
. To match the base R functionality, rmultinom_rcpp
returns a k
$\times$ n
IntegerMatrix
.
We now check if the rmultinom
and rmultinom_rcpp
functions produce the same results.
We generate a vector of 200 probabilities that sum to 1. We will sample 500 multinomial
outcomes and do this independently 20 times.
[1] TRUE
A benchmark of the functions suggests that the rmultinom_rcpp
function is very slightly
slower than the rmultinom
function, but that is not really a concern for our purposes.
Unit: milliseconds expr min lq mean median uq max neval cld rmultinom(1000, size, prob) 10.8676 10.9925 11.1737 11.0826 11.1898 13.9595 100 a rmultinom_rcpp(1000, size, prob) 11.0879 11.2072 11.4897 11.3341 11.6203 13.9920 100 b
The poisson distribution is a
non-negative discrete distribution characterized by having identical mean and
variance. The base R function rpois
used for generating Poisson data takes two
arguments: n
the number of simulated values to produce, and lambda
, a positive numeric
vector. The rpois
function cycles (and recycles) through the values in lambda
for each
successive value simulated. The function produces an integer vector of length n
. We
provide similar functionality using the R::rpois
function available in Rmath.h
. Note
that we cycle through the values of lambda
so that if the end of the lambda
vector is
reached before we have generated n
values, then we restart at the beginning of the
lambda
vector.
We now evaluate whether the rpois
and rpois
functions produce the same results. We
generate a positive vector with 200 values for lambda
and draw length(lambda) + 5
independent Poisson values.
[1] TRUE
A benchmark of the two functions suggests the rpois_rcpp
function may be slightly faster,
but once again, that is not our primary concern here.
Unit: microseconds expr min lq mean median uq max neval cld rpois(length(lambda) + 5, lambda) 7.412 7.7780 8.50529 7.9495 8.2165 60.014 100 b rpois_rcpp(length(lambda) + 5, lambda) 6.721 6.9965 7.37430 7.1645 7.4950 21.263 100 a
Thanks for visiting r-craft.org
This article is originally published at https://gallery.rcpp.org
Please visit source website for post related comments.