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