From 026342fd0e61d657312a9c75ebd03b3b9e59ed7c Mon Sep 17 00:00:00 2001 From: Bzero Date: Wed, 1 May 2024 09:04:57 +0200 Subject: [PATCH] Random distribution sampling (#416) * Add continuous uniform distribution sampling. * Add description for random(). * Add description for rand_uniform(). * Add discrete uniform distribution sampling * Add normal distribution sampling * Add geometric distribution sampling * Add exponential distribution sampling * Add log-normal distribution sampling * Add pareto distribution sampling * Add Bernoulli distribution sampling * Add Binomial distribution sampling * Add Poisson distribution sampling * Add random sampling functions to the list of predefined functions * Correct spelling Co-authored-by: David Peter * Turn rand_poisson() and rand_int() into scalar functions. * Use error() instead of NaN and make error messages more explicit in random sampling functions. --------- Co-authored-by: David Peter --- book/src/list-functions.md | 15 ++++++ numbat/modules/core/random.nbt | 99 ++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+) diff --git a/book/src/list-functions.md b/book/src/list-functions.md index 209d6c3e..61e090c8 100644 --- a/book/src/list-functions.md +++ b/book/src/list-functions.md @@ -101,6 +101,21 @@ fn sphere_area(radius: L) -> L^2 fn sphere_volume(radius: L) -> L^3 ``` +### Random sampling + +```nbt +fn random() -> Scalar +fn rand_uniform(a: T, b: T) -> T +fn rand_int(a: T, b: T) -> T +fn rand_normal(μ: T, σ: T) -> T +fn rand_bernoulli(p: Scalar) -> Scalar +fn rand_binomial(n: Scalar, p: Scalar) -> Scalar +fn rand_geometric(p: Scalar) -> Scalar +fn rand_poisson(λ: T) -> T +fn rand_exponential(λ: T) -> 1/T +fn rand_lognormal(μ: Scalar, σ: Scalar) -> Scalar +fn rand_pareto(α: Scalar, min: T) -> T + ### Algebra (experimental) ```nbt diff --git a/numbat/modules/core/random.nbt b/numbat/modules/core/random.nbt index ff81c3df..89ae6b8f 100644 --- a/numbat/modules/core/random.nbt +++ b/numbat/modules/core/random.nbt @@ -1,3 +1,102 @@ use core::scalar +use core::quantities +use core::error +use core::functions +use math::functions +# name: Standard uniform distribution sampling +# description: Uniformly samples the interval [0,1). fn random() -> Scalar + +# name: Continuous uniform distribution sampling +# url: https://en.wikipedia.org/wiki/Continuous_uniform_distribution +# description: Uniformly samples the interval [a,b) if a<=b or [b,a) if b(a: T, b: T) -> T = + if a <= b + then random() * (b - a) + a + else random() * (a - b) + b + +# name: Discrete uniform distribution sampling +# url: https://en.wikipedia.org/wiki/Discrete_uniform_distribution +# description: Uniformly samples the integers in the interval [a, b]. +fn rand_int(a: Scalar, b: Scalar) -> Scalar = + if a <= b + then floor( random() * (floor(b) - ceil(a) + 1) ) + ceil(a) + else floor( random() * (floor(a) - ceil(b) + 1) ) + ceil(b) + +# name: Bernoulli distribution sampling +# url: https://en.wikipedia.org/wiki/Bernoulli_distribution +# description: Samples a Bernoulli random variable, that is, 1 with probability p, 0 with probability 1-p. +# Parameter p must be a probability (0 <= p <= 1). +fn rand_bernoulli(p: Scalar) -> Scalar = + if p>=0 && p<=1 + then (if random() < p + then 1 + else 0) + else error("Argument p must be a probability (0 <= p <= 1).") + +# name: Binomial distribution sampling +# url: https://en.wikipedia.org/wiki/Binomial_distribution +# description: Samples a binomial distribution by doing n Bernoulli trials with probability p. +# Parameter n must be a positive integer, parameter p must be a probability (0 <= p <= 1). +fn rand_binom(n: Scalar, p: Scalar) -> Scalar = + if n >= 1 + then rand_binom(n-1, p) + rand_bernoulli(p) + else if n == 0 + then 0 + else error("Argument n must be a positive integer.") + +# name: Normal distribution sampling +# url: https://en.wikipedia.org/wiki/Normal_distribution +# description: Samples a normal distribution with mean μ and standard deviation σ using the Box-Muller transform. +fn rand_norm(μ: T, σ: T) -> T = + μ + sqrt(-2 σ² × ln(random())) × sin(2π × random()) + +# name: Geometric distribution sampling +# url: https://en.wikipedia.org/wiki/Geometric_distribution +# description: Samples a geometric distribution (the distribution of the number of Bernoulli trials with probability p needed to get one success) by inversion sampling. +# Parameter p must be a probability (0 <= p <= 1). +fn rand_geom(p: Scalar) -> Scalar = + if p>=0 && p<=1 + then ceil( ln(1-random()) / ln(1-p) ) + else error("Argument p must be a probability (0 <= p <= 1).") + +# A helper function for rand_poisson, counts how many samples of the standard uniform distribution need to be multiplied to fall below lim. +fn _poisson(lim: Scalar, prod: Scalar) -> Scalar = + if prod > lim + then _poisson(lim, prod × random()) + 1 + else -1 + +# name: Poisson distribution sampling +# url: https://en.wikipedia.org/wiki/Poisson_distribution +# description: Sampling a poisson distribution with rate λ, that is, the distribution of the number of events occurring in a fixed interval if these events occur with mean rate λ. +# The rate parameter λ must not be negative. +# This implementation is based on the exponential distribution of inter-arrival times. For details see L. Devroye, Non-Uniform Random Variate Generation, p. 504, Lemma 3.3. +fn rand_poisson(λ: Scalar) -> Scalar = + if λ >= 0 + then _poisson(exp(-λ), 1) + else error("Argument λ must not be negative.") + +# name: Exponential distribution sampling +# url: https://en.wikipedia.org/wiki/Exponential_distribution +# description: Sampling an exponential distribution (the distribution of the distance between events in a Poisson process with rate λ) using inversion sampling. +# The rate parameter λ must be positive. +fn rand_expon(λ: T) -> 1/T = + if value_of(λ) > 0 + then - ln(1-random()) / λ + else error("Argument λ must be positive.") + +# name: Log-normal distribution sampling +# url: https://en.wikipedia.org/wiki/Log-normal_distribution +# description: Sampling a log-normal distribution, that is, a distribution whose log is a normal distribution with mean μ and standard deviation σ. +fn rand_lognorm(μ: Scalar, σ: Scalar) -> Scalar = + exp( μ + σ × rand_norm(0, 1) ) + +# name: Pareto distribution sampling +# url: https://en.wikipedia.org/wiki/Pareto_distribution +# description: Sampling a Pareto distribution with minimum value min and shape parameter α using inversion sampling. +# Both parameters α and min must be positive. +fn rand_pareto(α: Scalar, min: T) -> T = + if value_of(min) > 0 && α > 0 + then min / ((1-random())^(1/α)) + else error("Both arguments α and min must be positive.")