# 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.