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, sample_method::DeterministicSamplingAlgorithm,
num_mats, T = Float64)
generate_design_matrices(n, d, sample_method::RandomSamplingAlgorithm,
num_mats, T = Float64)
generate_design_matrices(n, lb, ub, sample_method,
num_mats = 2)

Create num_mats matrices each containing a QMC point set, where:

n is the number of points to sample.
d is the dimensionality of the point set in [0, 1)ᵈ,
sample_method is the quasi-Monte Carlo sampling strategy used to create a deterministic point set out.
T is the eltype of the point sets. For some QMC methods (Faure, Sobol) this can be Rational

If the bound lb and ub are specified instead of d, the samples will be transformed into the box [lb, ub].

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, sample_method::DeterministicSamplingAlgorithm, num_mats, T = Float64)

Create an iterator for doing multiple i.i.d. randomization over QMC sequences where

num_mats is the lenght of the iterator
n is the number of points to sample.
d is the dimensionality of the point set in [0, 1)ᵈ,
sample_method is the quasi-Monte Carlo sampling strategy used to create a deterministic point set out.
T is the eltype of the point sets. For some QMC methods (Faure, Sobol) this can be Rational

It is now compatible with all scrambling methods and shifting. One can also use it with Distributions.Sampleable or RandomSample.

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.