Distributions

Our distribution support will grow over time. If you need a distribution that is not yet supported please post an issue <https://github.com/mrc-ide/dust/issues/>.

Uniform

template<typename real_type, typename rng_state_type>
__nv_exec_check_disable__ __host__ __device__ real_type dust::random::uniform(rng_state_type &rng_state, real_type min, real_type max)

Draw a uniformly distributed random number with arbitrary bounds. This function simply scales the output of dust::random::random_real

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • min, max – The minimum and maximum values for the distribution

Normal

template<typename real_type, algorithm::normal algorithm = algorithm::normal::box_muller, typename rng_state_type>
__nv_exec_check_disable__ __host__ __device__ real_type dust::random::normal(rng_state_type &rng_state, real_type mean, real_type sd)

Draw a normally distributed random number with arbitrary bounds. This function simply scales the output of dust::random::random_real

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • The – algorithm to use; the default is Box-Muller which is slowest on CPU but simple.

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • mean – The mean of the distribution

  • sd – The standard deviation of the distribution

The enum dust::random::algorithm::normal controls the normal algorithm used for a given draw

enum class dust::random::algorithm::normal

Values:

enumerator box_muller

Box-Muller method (transformation)

enumerator polar

Polar method (rejection)

enumerator ziggurat

Ziggurat method (rejection)

Binomial

template<typename real_type, typename rng_state_type>
__host__ __device__ real_type dust::random::binomial(rng_state_type &rng_state, real_type n, real_type p)

Draw a binomially distributed random number; the number of successes given n trials each with probability p. Generation is performed using a rejection-sampling algorithm or inversion depending on the expected mean (n * p).

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • n – The number of trials

  • p – The probability of success of each trial

Hypergeometric

template<typename real_type, typename rng_state_type>
__host__ __device__ real_type dust::random::hypergeometric(rng_state_type &rng_state, real_type n1, real_type n2, real_type k)

Draw random number from the hypergeometric distribution. This is often descrribed as sampling k elements without replacement from an urn with n1 white balls and n2 black balls, where the outcome is the number of white balls retrieved.

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • n1 – The number of white balls in the urn

  • n2 – The number of black balls in the urn

  • k – The number of draws

Exponential

template<typename real_type, typename rng_state_type>
__nv_exec_check_disable__ __host__ __device__ real_type dust::random::exponential(rng_state_type &rng_state, real_type rate)

Draw a exponentially distributed random number given a rate parameter. Generation is performed using inversion (faster algorithms exist but are not yet implemented).

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • rate – The rate of the process

Poisson

template<typename real_type, typename rng_state_type>
__nv_exec_check_disable__ __host__ __device__ real_type dust::random::poisson(rng_state_type &rng_state, real_type lambda)

Draw a Poisson distributed random number given a mean parameter. Generation is performed using either Knuth’s algorithm (small lambda) or Hormann’s rejection sampling algorithm (medium lambda), or rejection based on the Cauchy (large lambda)

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • lambda – The mean of the distribution

Returns:

We return a real_type, not an int, as with deterministic mode this will not necessarily be an integer

Multinomial

template<typename real_type, typename rng_state_type, typename T, typename U>
__host__ __device__ void dust::random::multinomial(rng_state_type &rng_state, int size, const T &prob, int prob_len, U &ret)

Draw one sample from the multinomial distribution.

This is written assuming that prob and ret are arbitrary containers; they could be pointers, vectors or anything else that supports random access. We could also do this with iterators but they always feel a bit weird really.

Note that prob and ret will ordinarily be containers of real_type but that might not be the case, in particular when calling from R where we want calculations to happen in single precision but the inputs (and destination) are double precision.

Later we will provide some helpers that will allow setting a stride over these containers; this will help with use from odin code.

Template Parameters:
  • real_type – The underlying real number type, typically double or float. A compile-time error will be thrown if you attempt to use a non-floating point type (based on `std::is_floating_point).

  • rng_state_type – The random number state type

  • T, U – The type of the containers for prob and ret. This might be double* or std::vector<double> depending on use.

Parameters:
  • rng_state – Reference to the random number state, will be modified as a side-effect

  • size – The number of trials (analagous to size in binomial)

  • prob – The set of probabilities. In a binomial trial we only provide a single probability of success but here there is a vector/array of such probabilities. These need not sum to one as they will be normalised, however they must all be non-negative and at least one must be positive.

  • prob_len – The number of probabilities (or outcomes)

  • ret – Container for the return value

Example

#include <iostream>
#include <dust/random/random.hpp>

void show_numbers(std::vector<double> numbers, const char * label) {
  std::cout << label << ":";
  for (size_t i = 0; i < numbers.size(); ++i) {
    std::cout << " " << numbers[i];
  }
  std::cout << std::endl;
}

int main() {
  using rng_state_type = dust::random::generator<double>;

  auto state = dust::random::seed<rng_state_type>(42);
  size_t n = 10;

  // A place to store numbers:
  std::vector<double> numbers(n);

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::uniform<double>(state, 10, 20);
  }
  show_numbers(numbers, "Uniform(10, 20)");

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::normal<double>(state, 5, 2);
  }
  show_numbers(numbers, "Normal(5, 2)");

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::exponential<double>(state, 6);
  }
  show_numbers(numbers, "Exponential(6)");

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::poisson<double>(state, 3);
  }
  show_numbers(numbers, "Poisson(3)");

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::binomial<double>(state, 10, 0.3);
  }
  show_numbers(numbers, "Binomial(10, 0.3)");

  for (size_t i = 0; i < n; ++i) {
    numbers[i] = dust::random::hypergeometric<double>(state, 10, 3, 7);
  }
  show_numbers(numbers, "Hypergeometric(10, 3, 7)");
}
Uniform(10, 20): 14.9695 16.8962 13.6696 13.6204 19.2999 10.7353 13.6703 10.8023 16.4116 18.2632
Normal(5, 2): 5.02109 6.01904 5.15919 6.38409 2.67724 6.1001 4.48944 5.31519 3.86882 5.23517
Exponential(6): 0.0231662 0.0920166 0.238532 0.0190538 0.241245 0.128207 0.13686 0.0170973 0.0575175 0.097848
Poisson(3): 2 0 2 2 4 3 4 3 2 4
Binomial(10, 0.3): 3 4 3 1 4 2 3 4 4 3
Hypergeometric(10, 3, 7): 6 6 6 6 7 5 6 6 5 6