Distributions

Loading Distributions activates fixed-work methods for a small distribution set.

using PureRNGs, Random, Distributions

rng = Philox4x32(123456)
distribution = LogNormal(2.0, 0.5)
values, rng = rand_next(rng, distribution, 1024)

buffer = similar(values)
_, rng = rand_next!(rng, distribution, buffer)
x = rand_at(rng, distribution, 1)

Direct immutable methods

Supported fixed-work distributions are Normal, Uniform, Exponential, LogNormal, Weibull, Rayleigh, Laplace, Logistic, Gumbel, Pareto, Frechet, Cauchy, TriangularDist, Bernoulli, and DiscreteUniform.

The continuous distributions and Bernoulli require Float32 or Float64 parameters. Bernoulli returns Bool. DiscreteUniform returns Int. Other destinations must match the distribution's result type.

Each distribution supports scalar, addressed, allocating, and destination-fill draws. Every sampling form also has a continuation form where applicable.

Parameters must be finite. Normal and LogNormal permit zero scale. Uniform requires increasing bounds with a finite difference. Exponential and Rayleigh require positive scale. Weibull requires positive shape and scale. Laplace requires a positive scale. Bernoulli requires a probability in [0, 1]. Logistic, Gumbel, and Cauchy require positive scales. Pareto and Frechet require positive shapes and scales. TriangularDist requires a ≤ c ≤ b and a finite width b-a; endpoint modes and a == b == c are supported.

The fixed transforms and their spans per result are:

Here u uses the same open midpoint grid as normal draws, before the normal transform. It never equals zero or one. Gumbel and Frechet use e = -log(1-u) with the generator's exponential backend transform, or native compiled math under Reactant. This avoids singular input endpoints without retries. TriangularDist uses an ordinary uniform draw and a normalized inverse CDF that avoids products of interval widths. Even a point mass consumes its full span.

Exponential draws are strictly positive: the smallest is 5.96e-8 for Float32 and 1.11e-16 for Float64. For comparison, Base's ziggurat randn reaches 13.708 and randexp reaches 44.434.

Float32 normals stop at 5.29 sigma and drop one draw in 8.4 million, so tail-sensitive Float32 work should draw Float64 and convert.

Every derived distribution inherits these bounds through its transform. Bernoulli(p) realizes ceil(p * 2^w) / 2^w with w equal to 24 or 53, so probabilities below 2^-w round up to 2^-w. Use Float64 parameters for p below about 1e-7. Categorical shares the weighted sampling resolution described in Sampling.

distribution = Logistic(0.0, 1.0)
values, rng = rand_next(rng, distribution, 128)
_, rng = rand_next!(rng, distribution, values)

LogNormal(μ, 0) still consumes its normal span. Valid transforms may underflow or overflow under IEEE arithmetic; PureRNGs neither clips values nor retries a draw. The exact random bits, order, and continuation agree within an execution site, while transcendental transforms may differ between CPU, CUDA, and compiled execution.

CUDA supports these direct methods. AMDGPU serves the same direct methods as CUDA. Reactant supports scalar, continuation, and addressed forms for these continuous distributions, using its native compiled arithmetic. This does not add parameter-gradient support or a general distribution AD guarantee. Metal serves every distribution with Float32 parameters; Float64 parameters need Float64 arithmetic, which Metal lacks.

Categorical

Categorical draws return Int labels in 1:length(p). They use the raw probability values as weights over those labels and the same 53-bit Float64 weighted mapping as randsample; there is no extra normalization or alias table. The probabilities must convert to finite, nonnegative Float64 values with a finite, positive cumulative total. CUDA's parallel preparation can change rounding and selected categories compared with the CPU, as described in Supply weights.

rng = Philox4x32(123456)
distribution = Categorical([0.1, 0.7, 0.2])
labels, rng = rand_next(rng, distribution, 1024)

buffer = similar(labels)
_, rng = rand_next!(rng, distribution, buffer)
label = rand_at(rng, distribution, 1)

Scalar and addressed Categorical calls validate their weights on the CPU and accept only CPU-backed or device-agnostic probabilities. Small populations use an allocation-free scan; larger populations build a cumulative table. Both observe changed probabilities on the next call. Batch and destination-fill calls prepare once per call: CPU accepts CPU-backed or device-agnostic probabilities, while a GPU requires probabilities on the matching device. A mismatch throws ArgumentError rather than copying data. Categorical is not a Reactant carrier operation and adds no AD support. Metal rejects device-executing allocation and fills, since the cumulative table is Float64.

MvNormal

Fill destinations must not share storage with the mean or covariance factor.

An MvNormal draw applies Distributions' own map, μ + L z, to the next length(d) standard normal draws, where L is the covariance's lower factor. n draws are one normal fill, a column per draw, so a matrix of draws equals the chained single draws.

rng = Philox4x32(123456)
distribution = MvNormal([1.0, -2.0], [2.0 0.3; 0.3 1.0])
point, rng = rand_next(rng, distribution)
cloud, rng = rand_next(rng, distribution, 1000)
seventh = rand_at(rng, distribution, 7)

PDMat, diagonal, and scalar covariances work on every backend. Sparse PDSparseMat covariances support allocating and in-place CPU draws. A device draw moves the covariance's factor to the device once per call and applies it there. Gradients with respect to μ and the covariance are pathwise with Enzyme, Mooncake, and ForwardDiff.

Gamma and its family

Gamma, Chisq, InverseGamma, Beta, TDist, and Dirichlet draws use the rejection sampler of Marsaglia and Tsang (2000) with a fixed stream span. Each Gamma draw reserves one boost uniform and eight candidates, each a normal and an open uniform, and takes the first candidate the test accepts. If all eight reject, which happens about once in 3·10¹⁰ draws at shape 1 and far less often at larger shapes, the draw continues the same test on a child stream keyed by its position. The fallback preserves the rejection rule, and the parent stream advances by the same 17 normal widths either way. Shapes below one draw Gamma(shape + 1) and multiply by u^(1/shape) from the boost uniform.

rng = Philox4x32(123456)
waits, rng = rand_next(rng, Gamma(2.5, 2.0), 1024)
proportion = rand_at(rng, Beta(0.3, 0.4), 7)
mixture, rng = rand_next(rng, Dirichlet([0.3, 1.0, 2.5]))

Chisq(ν) is Gamma(ν/2, 2), and InverseGamma(α, θ) is θ over a Gamma(α) draw. When a tiny shape makes the standard draw subnormal or zero, a log-space repair preserves representable scaled results. It rereads the same bits and leaves parent advancement unchanged. Beta and Dirichlet normalize log-Gamma draws in consecutive spans. When tiny positive shapes overflow those logs, a scaled comparison reuses the same bits to avoid subtracting infinities. Results can round to zero or one. Dirichlet columns remain normalized, and the parent span stays unchanged.

A Dirichlet fill destination must not share storage with its concentration vector. Host concentration arrays transfer to the generator's backend when needed. Concentrations already on that backend stay there, including device views. Concentrations on another backend are rejected rather than copied through the CPU. TDist(ν) divides a normal by the square root of a following chi-square over ν. Every member runs on the CPU, CUDA, AMDGPU, and, with Float32 parameters, Metal; a 2^26 Gamma fill takes about 6.7 ms on an A100. CUDA Dirichlet fills use one work item per column when enough columns fill the kernel's resident blocks. Smaller column batches and the other device backends fill components in parallel, then normalize each column in component order. Under Reactant the scalar members have scalar, continuation, and addressed forms: a compiled draw evaluates all eight candidates, selects the first accepted one without branching, and runs the child stream as a traced loop.

Shape gradients use the implicit derivative of the Gamma CDF, dx/dα = -∂F(x; α)/∂α ÷ f(x; α) (Figurnov, Mohamed, and Mnih 2018), with Enzyme, Mooncake, and ForwardDiff alike. Differentiating through the rejection test would bias the gradient, so the sampler is a primitive with that rule. The derivative uses differentiated incomplete-gamma series and continued fractions. Large shapes near the mean use a constant-work central expansion, avoiding cancellation and stalled floating-point increments. The arithmetic runs in the shape's type, including inside device kernels. The central branch begins at shape 10⁴ for Float32 and 10⁶ for Float64. Finite draws do not guarantee finite derivatives at extreme shapes: derivative arithmetic and upstream distribution constructors can overflow. Rounded endpoints have zero pathwise tangents where the sampler controls differentiation. In particular, Mooncake can produce NaN gradients for subnormal concentrations even when the sampled values are finite.

Other distributions

Use the mutable bridge for Distributions' wider API:

rng = StatefulRNG(Philox4x32(123456))
values = rand(rng, Poisson(3.0), 1024)

This uses Distributions' algorithms, not the package's fixed-work mappings. Even supported distribution names can produce different streams through the bridge.