Sampler APIs

sample(n::Integer, d::Integer, sampler::SamplingAlgorithm, T = Float64)
sample(n::Integer, lb, ub, sampler::SamplingAlgorithm)

Generate n sample points with sampler.

Arguments

• n: Positive number of points to generate.
• d: Positive dimension of the unit box. This form returns a d-by-n matrix with values in [0, 1].
• lb: Scalar, tuple, or vector of lower bounds. Its length determines the dimension when it is not scalar.
• ub: Scalar, tuple, or vector of upper bounds with the same shape as lb. Each lower bound must be less than or equal to its corresponding upper bound.
• sampler: Concrete SamplingAlgorithm that determines the point set construction.

Arguments with defaults

• T = Float64: Element type of the unit-box result. The bounds form infers its output element type from lb.

Returns

A matrix whose columns are the generated points. The unit-box form has size (d, n). The bounds form has size (length(lb), n) for collection bounds and maps every coordinate from [0, 1] to its corresponding closed interval [lb[i], ub[i]].

Throws

• AssertionError: If n is not positive, the bounds have different lengths, or a lower bound exceeds its corresponding upper bound.

Examples

julia> using QuasiMonteCarlo

julia> unit_points = QuasiMonteCarlo.sample(4, 2, SobolSample());

julia> size(unit_points)
(2, 4)

julia> points = QuasiMonteCarlo.sample(4, [0.0, -1.0], [1.0, 1.0], SobolSample());

julia> all([0.0, -1.0] .<= points .<= [1.0, 1.0])
true

Extension rules

To add a sampler, subtype SamplingAlgorithm and implement only the unit-box form sample(n, d, sampler, T). The implementation must return a d-by-n matrix with elements in [0, 1]; the bounds form is provided by this package and delegates to that unit-box method. Do not extend this bounds method for a new sampler.

SamplingAlgorithm

Abstract supertype for sampling strategies accepted by sample.

Extension rules

Define a concrete subtype and implement sample(n::Integer, d::Integer, sampler::YourSampler, T = Float64). The method must return a d-by-n matrix of elements of type T whose columns are points in the unit box [0, 1]^d. It must reject invalid sampler-specific inputs with an informative exception and must preserve the requested number of points and dimensions. sample applies user-supplied bounds itself, so extensions should implement the unit-box method rather than a separate bounds method.

GridSample(R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

A simple rectangular grid lattice. It samples from a regular grid in the unit box and optionally applies R.

Fields

• R::RandomizationMethod = NoRand(): Randomization applied to the grid.

In more than 2 dimensions, grids have worse discrepancy than simple Monte Carlo. As a result, they should almost never be used for multivariate integration; their use is as a starting point for other algorithms.

SobolSample(R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

Samples taken from Sobol's base-2 sequence.

Fields

• R::RandomizationMethod = NoRand(): Randomization applied to the Sobol points.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, SobolSample());

julia> size(points)
(2, 8)

VanDerCorputSample(base::Integer, R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

The van der Corput sequence, also called the radical inverse sequence, is a one-dimensional low-discrepancy sequence that recursively splits the unit interval into equally-sized pieces before inserting one sample in each interval.

For example, in base 2, the sequence starts by inserting one point in each half of the unit interval in the first pass (the first 2 samples); then one sample in each quarter, then each eighth, and so on. This creates a well-stratified sample, so long as the number of samples is a multiple of a power of the base. The one-dimensional result is returned as a vector of length n.

Fields

• base::Integer: Base used for the radical-inverse expansion.
• R::RandomizationMethod = NoRand(): Randomization applied to the sequence.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 1, VanDerCorputSample(base = 2));

julia> size(points)
(8,)

FaureSample(R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

A Faure low-discrepancy sequence in the prime base determined by the sample dimension.

Fields

• R::RandomizationMethod = NoRand(): Randomization applied after the deterministic Faure points are generated.

Faure-distributed samples cover all dimensions evenly, using the same set of points for all variables, up to ordering.

When scrambled, randomized Faure sequences provide worst-case guarantees that variance will be at most exp(1) ≈ 2.718 times greater than for a purely Monte Carlo integral. However, they are much less efficient than the Sobol sequence at integrating functions with low effective dimension (functions where the first few inputs dominate the evaluation).

The Faure sequence in dimension s forms a (0, s)-sequence with base b = nextprime(s).

A Faure sequence must have a supported length of k * base^s with an integer k satisfying 1 ≤ k < base.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, FaureSample());

julia> size(points)
(2, 8)

References: Faure, H. (1982). Discrépance de suites associées à un système de numération (en dimension s). Acta Arith., 41, 337-351. Owen, A. B. (1997). Monte Carlo variance of scrambled net quadrature. SIAM Journal on Numerical Analysis, 34(5), 1884-1910.

LatticeRuleSample(R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

Generate a point set using a rank-1 lattice rule.

Fields

• R::RandomizationMethod = NoRand(): Randomization applied to the lattice points.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, LatticeRuleSample());

julia> size(points)
(2, 8)

HaltonSample(R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm

Create a Halton sequence using successive prime bases.

Fields

• R::RandomizationMethod = NoRand(): Randomization applied to the sequence.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, HaltonSample());

julia> size(points)
(2, 8)

GoldenSample()

Generate a quasirandom Kronecker sequence using powers of the generalized golden ratio.

The harmonious, or generalized golden, ratios are defined as the solutions to the equation: $x^d = x + 1$

Where d is the dimension of the sequence. The Golden sequence is then equivalent to Kronecker([x^-i for i in 1:d]).

WARNING: the generalized golden sequence in more than 2 dimensions is purely experimental. It likely has poor discrepancy in high dimensions, and should not be used without verifying answers against a better-known quasirandom sequence. Try a rank-1 lattice rule instead.

References: Roberts, M. (2018). The Unreasonable Effectiveness of Quasirandom Sequences. Extreme Learning. http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/

Returns

• A KroneckerSample configured with a generalized golden-ratio generator.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, GoldenSample());

julia> size(points)
(2, 8)

KroneckerSample(generator::AbstractVector, R::RandomizationMethod = NoRand()) <: DeterministicSamplingAlgorithm
KroneckerSample(d::Integer, R::RandomizationMethod = NoRand(), T = Float64)

A Kronecker sequence is a point set generated using a vector and the equation x[i] = i * generator .% 1

Where i runs from 1 through the sample size n. This sequence will be equidistributed (uniform in the infinite limit) so long as the components of generator are linearly independent over the field of rational numbers.

If no generator is specified, a lattice based on the generalized golden ratio is used; see GoldenSample for more information.

Kronecker sequences are not recommended for use in more than 3 dimensions, as theory on them is sparse. LatticeRuleSample will return rank-1 lattice rules, which behave similarly to Kronecker sequences but have better properties.

Fields

• generator::AbstractVector: Generator vector. Its length must equal the requested dimension.
• R::RandomizationMethod = NoRand(): Randomization applied to the generated sequence.

Arguments

• d::Integer: Dimension for the constructor that creates a generalized golden-ratio generator.
• R::RandomizationMethod: Optional randomization method.
• T: Element type used for the generated generator.

Examples

julia> using QuasiMonteCarlo

julia> points = sample(8, 2, KroneckerSample(2));

julia> size(points)
(2, 8)

References: Leobacher, G., & Pillichshammer, F. (2014). Introduction to quasi-Monte Carlo integration and applications. Switzerland: Springer International Publishing. https://link.springer.com/content/pdf/10.1007/978-3-319-03425-6.pdf

RandomSample(; rng = Random.TaskLocalRNG()) <: RandomSamplingAlgorithm

Plain Monte Carlo sampler that draws independent uniform samples from the unit box.

Fields

• rng::AbstractRNG: random-number generator used for each call to sample.

Keywords

• rng::AbstractRNG = Random.TaskLocalRNG(): random-number generator used to draw the samples. Supply a seeded RNG for reproducible results.

Examples

julia> using QuasiMonteCarlo, Random

julia> sampler = RandomSample(MersenneTwister(42));

julia> points = sample(4, 2, sampler);

julia> size(points)
(2, 4)

LatinHypercubeSample(rng::AbstractRNG = Random.TaskLocalRNG()) <: RandomSamplingAlgorithm

A Latin hypercube is a point set with the property that every one-dimensional interval (i / n, (i + 1) / n) contains exactly one point. It is useful for high-dimensional sampling because it is more uniform than independent Monte Carlo sampling without requiring a full tensor-product grid.

Fields

• rng::AbstractRNG = Random.TaskLocalRNG(): Random-number generator used to independently permute the strata in each dimension.

Examples

julia> using QuasiMonteCarlo, Random

julia> points = sample(8, 2, LatinHypercubeSample(MersenneTwister(42)));

julia> size(points)
(2, 8)

RandomizedHaltonSample(; rng = Random.TaskLocalRNG()) <: RandomSamplingAlgorithm

Create a randomized Halton sequence by independently permuting the radical inverse digits in each dimension.

Fields

• rng::AbstractRNG = Random.TaskLocalRNG(): Random-number generator used for the digit permutations.

Examples

julia> using QuasiMonteCarlo, Random

julia> points = sample(8, 2, RandomizedHaltonSample(rng = MersenneTwister(42)));

julia> size(points)
(2, 8)

References: Owen, A. (2017). A randomized Halton algorithm in R. https://doi.org/10.48550/arXiv.1706.02808