Create new generators

Summary

This document describes how you can extend RcppRNG with your own generator classes for random numbers.

As an example, we will create an implementation of a generator for the Erlang distribution, see https://en.wikipedia.org/wiki/Erlang_distribution. Let \(k \in \mathbb{N}\) and \(\lambda > 0\). The \(\mathrm{Erlang}{(k, \lambda)}\) distribution has the stochastic representation \[ \sum_{j=1}^{k} E_j , \] where \(E_1, \ldots, E_k\) are iid exponentially distributed with rate \(\lambda\).

Setup

First, you have to tell Rcpp that you are using c++11 and that your code links to the packages RcppRNG (for obvious reasons) as well as dqrng, BH, and sitmo (to use the RNG from dqrng).

#include <Rcpp.h>
// [[Rcpp::plugins("cpp14")]]
// [[Rcpp::depends(RcppRNG,dqrng,BH,sitmo)]]

#include <RcppRNG.hpp>

// for dqrng functions
#include <R_randgen.h>
#include <convert_seed.h>
#include <dqrng_generator.h>
#include <dqrng_distribution.h>
#include <xoshiro.h>

Implementation

The distribution

namespace RcppRNG {

namespace distribution {

namespace erlang {

template <typename _IntType, typename _RealType>
void check_params(const _IntType shape, const _RealType lambda) {
  if (!(0 < shape)) throw std::range_error("shape not > 0");
  if (!(0. <= lambda && lambda <= std::numeric_limits<_RealType>::infinity()))
    throw std::range_error("lambda not in [0, infty]");
}

}  // namespace erlang

template <typename _IntType, typename _RealType,
          typename _ExponentialDistributionType>
class erlang_distribution {
 public:
  using result_type = _RealType;
  static_assert(std::is_integral<_IntType>::value,
                "_IntType not integral type");
  static_assert(std::is_floating_point<_RealType>::value,
                "_RealType not floating point number");
  static_assert(std::numeric_limits<_RealType>::is_iec559,
                "_RealType not IEEE 754");

  class param_type {
   public:

    explicit param_type(const _IntType shape = _IntType{1},
                        const _RealType lambda = _RealType{1.})
        : shape_{shape}, lambda_{lambda} {
      erlang::check_params(shape_, lambda_);
    }

    // Used for construction from a different specialization
    template <typename _ErlangParamType,
              class = typename std::enable_if<
                  !std::is_convertible<_ErlangParamType, param_type>::value>>
    explicit param_type(const _ErlangParamType& params)
        : shape_{params.shape()}, lambda_{params.lambda()} {}

    // compiler generated ctor and assignment op is sufficient

    _IntType shape() const { return shape_; }

    _RealType lambda() const { return lambda_; }

   private:
    _IntType shape_;
    _RealType lambda_;
  };

  explicit erlang_distribution(const _IntType shape = _IntType{1},
                               const _RealType lambda = _RealType{1.})
      : shape_{shape}, exponential_distribution_{lambda} {
    erlang::check_params(shape_, exponential_distribution_.lambda());
  }

  explicit erlang_distribution(const param_type& params)
      : shape_{params.shape()}, exponential_distribution_{params.lambda()} {}

  // compiler generated ctor and assignment op is sufficient

  param_type params() const {
    return param_type{shape_, exponential_distribution_.lambda()};
  }

  void params(const param_type& params) {
    shape_ = params.shape();
    exponential_distribution_.params(params.lambda());
  }

  _IntType shape() const { return shape_; }

  _RealType lambda() const { return exponential_distribution_.lambda(); }

  template <typename _EngineType>
  result_type operator()(_EngineType&& engine) {
    result_type out = _RealType{0.};
    for (_IntType j = _IntType(1); j <= shape_; ++j)
      out += exponential_distribution_(std::forward<_EngineType>(engine));

    return out;
  }

  template <typename _EngineType>
  result_type operator()(_EngineType&& engine, const param_type& params) {
    return typename param_type::distribution_type{params}(
        std::forward<_EngineType>(engine));
  }

 private:
  _IntType shape_;
  _ExponentialDistributionType exponential_distribution_;
};  // erlang_distribution

}  // namespace distribution

template <typename _IntType, typename _RealType,
          typename _ExponentialDistributionType>
using erlang_distribution =
    distribution::erlang_distribution<_IntType, _RealType,
                                      _ExponentialDistributionType>;

}  // namespace RcppRNG

Rcpp-wrapper

using namespace Rcpp;

using unit_exponential_distribution =
    RcppRNG::unit_exponential_distribution<double>;
using exponential_distribution =
    RcppRNG::exponential_distribution<double, unit_exponential_distribution>;
using erlang_distribution =
    RcppRNG::erlang_distribution<int, double, exponential_distribution>;

//' @export
// [[Rcpp::export(rng=false)]]
NumericVector Rcpp_rerlang_RcppRNG(R_xlen_t n, int shape = 1, double rate = 1.) {
  NumericVector out{no_init(n)};
  RcppRNG::rng::RcppRNG engine{};
  erlang_distribution erlang{shape, rate};
  std::generate(out.begin(), out.end(),
                [&engine, &erlang]() { return erlang(engine); });

  return out;
}

//' @export
// [[Rcpp::export(rng=false)]]
NumericVector Rcpp_rerlang_DQRNG(R_xlen_t n, int shape = 1, double rate = 1.) {
  NumericVector out{no_init(n)};
  RcppRNG::rng::DQRNG engine{};
  erlang_distribution erlang{shape, rate};
  std::generate(out.begin(), out.end(),
                [&engine, &erlang]() { return erlang(engine); });

  return out;
}

Test

set.seed(1623L)
Rcpp_rerlang_RcppRNG(10, 3, 0.5)
#>  [1] 10.502098  7.088853  5.084991  7.435394  2.082935  7.849784  5.739142
#>  [8]  9.858097  8.521792  6.443745
RcppRNG::dqset_seed(1623L)
Rcpp_rerlang_DQRNG(10, 3, 0.5)
#>  [1]  3.401419  5.063018  4.915893 11.271596 14.128176  7.942866  8.195575
#>  [8]  2.825962 12.446754  3.054656

Benchmark

n <- 1e5L
shape <- 10L
lambda <- 0.5
bm <- bench::mark(
  Rcpp = Rcpp_rerlang_RcppRNG(n, shape, lambda),
  dqrng = Rcpp_rerlang_DQRNG(n, shape, lambda),
  check=FALSE,
  min_iterations = 20L
)
ggplot2::autoplot(bm)
#> Loading required namespace: tidyr