Design Matrices

API

It is often convenient to generate multiple independent sequences, for error estimation (uncertainty quantification). The resulting sequences can be stored in what is often called a design matrix. In this package, this is achieved with the generate_design_matrices(n, d, ::DeterministicSamplingAlgorithm), ::RandomizationMethod, num_mats) function. num_mats is the number of independent realizations. The resulting design matrix is a vector of matrix of length num_mats.

QuasiMonteCarlo.generate_design_matrices — Function

generate_design_matrices(n, d, sampler, num_mats, T = Float64)
generate_design_matrices(n, lb, ub, sampler, num_mats = 2)

Generate multiple point-set matrices with sampler.

Arguments

n: Positive number of points in each generated matrix.
d: Positive dimension of the unit box. This form returns matrices of size (d, n) with elements in [0, 1].
lb: Collection of lower bounds. Its length determines the dimension.
ub: Collection of upper bounds with the same length as lb.
sampler: Concrete SamplingAlgorithm used to construct each point set.
num_mats: Number of matrices to generate.

Optional positional arguments

T = Float64: Element type of the unit-box matrices.

Returns

A vector of num_mats matrices. Each unit-box matrix has size (d, n). The bounds form maps every coordinate to its corresponding interval [lb[i], ub[i]].

Throws

AssertionError: If the bounds have different lengths or a lower bound exceeds its corresponding upper bound.
An exception from the selected sampler if n, d, or num_mats is invalid.

Examples

julia> using QuasiMonteCarlo

julia> matrices = generate_design_matrices(4, 2, RandomSample(), 3);

julia> length(matrices)
3

julia> all(size(matrix) == (2, 4) for matrix in matrices)
true

Developer Notes

This function builds on the public sample contract. New sampling algorithms should extend sample; they should not add methods to generate_design_matrices.

source

generate_design_matrices(n, d, sampler, R::NoRand, num_mats, T = Float64)

R = NoRand() produces num_mats matrices each containing a different deterministic point set in [0, 1)ᵈ. Note that this is an ad hoc way to produce i.i.d sequence as it creates a deterministic point in dimension d × num_mats and split it in num_mats point set of dimension d. This does not have any QMC garantuees.

source

Instead of generating num_mats matrices, it is possible (and more memory efficient) to randomize the same matrix multiple times and perform an operation after each randomization. This is possible using iterators. To build an iterator use the DesignMatrix function.

QuasiMonteCarlo.DesignMatrix — Function

DesignMatrix(n, d, sampler, num_mats, T = Float64)
DesignMatrix(n, d, sampler, randomization, num_mats, T = Float64)

Create an iterator that produces num_mats point-set matrices.

The iterator is useful when a computation needs independent randomized QMC realizations but allocating all matrices at once is undesirable. For a deterministic sampler, the five-argument form uses the sampler's R field; the six-argument form selects randomization explicitly. RandomSample also has a five-argument form and generates independent Monte Carlo matrices.

Arguments

n::Integer: Number of points in each matrix.
d::Integer: Dimension of each point. Each yielded matrix has size (d, n).
sampler: Deterministic or random SamplingAlgorithm.
randomization::RandomizationMethod: Randomization applied to a deterministic sampler.
num_mats::Integer: Number of matrices yielded by the iterator.
T::DataType = Float64: Element type of each matrix.

Returns

An AbstractDesignMatrix iterator whose elements are matrices of size (d, n) and element type T.

Throws

An exception from the selected sampler or randomization method if n, d, the randomization base, or the sampler-specific parameters are invalid.

Examples

julia> using QuasiMonteCarlo

julia> iterator = DesignMatrix(4, 2, SobolSample(R = Shift()), 3);

julia> length(iterator)
3

julia> size(first(iterator))
(2, 4)

Extension rules

For a custom deterministic sampler, define the unit-box sample method from SamplingAlgorithm and provide an R::RandomizationMethod field. The built-in randomization-specific initialize methods are implementation details; custom extensions should use the public sample and randomize contracts unless they are deliberately implementing a compatible iterator.

source

Warning

The method generate_design_matrices(n, d, sampler, R::NoRand, num_mats, T = Float64) is an ad hoc way to produce a Design Matrix. Indeed, it creates a deterministic point set in dimension d × num_mats and splits it into num_mats point set of dimension d. The resulting sequences have no QMC guarantees. This seems to have been proposed in Section 5.1 of Saltelli, A. et al. (2010)^[1] to do uncertainty quantification. See this discussion for a visual proof.

Example

using QuasiMonteCarlo, Random, StatsBase
Random.seed!(1234)
m = 4
d = 3
b = QuasiMonteCarlo.nextprime(d)
N = b^m # Number of points
pad = 2m # Can also choose something as `2m` to get "better" randomization
num_mats = 5

f(x) = prod(x) * 2^length(x) # test function ∫f(x)dᵈx = 1

# Randomize over num_mats = 5 independent Randomized Faure sequences
iterator = DesignMatrix(N, d, FaureSample(R = OwenScramble(base = b, pad = pad)), num_mats)

μ = [mean(f(c) for c in eachcol(X)) for X in iterator]

5-element Vector{Float64}:
 0.993319147196098
 0.997134058324563
 1.0054957094324526
 0.9923122929309586
 1.001184984852591

Using std(μ) then gives you the estimated variance of your RQMC prediction.

# Or using `generate_design_matrices`. Note that this is less memory efficient since it allocate space for 5 large big matrices.
μ = [mean(f(c) for c in eachcol(X))
     for X in QuasiMonteCarlo.generate_design_matrices(N,
    d,
    FaureSample(R = OwenScramble(base = b, pad = pad)),
    num_mats)]

5-element Vector{Float64}:
 0.9848852316008134
 0.9909436902184221
 1.0108023981060077
 1.009804779326848
 0.9831653609748305

1Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., & Tarantola, S. (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer physics communications, 181(2), 259-270.