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
orfloat
. 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
orfloat
. 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
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 probabilityp
. 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
orfloat
. 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 withn1
white balls andn2
black balls, where the outcome is the number of white balls retrieved.- Template Parameters:
real_type – The underlying real number type, typically
double
orfloat
. 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
orfloat
. 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
orfloat
. 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 anint
, 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
andret
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
andret
will ordinarily be containers ofreal_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
orfloat
. 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
andret
. This might bedouble*
orstd::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
inbinomial
)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