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
— Functiongenerate_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 setout
.T
is theeltype
of the point sets. For some QMC methods (Faure, Sobol) this can beRational
If the bound lb
and ub
are specified instead of d
, the samples will be transformed into the box [lb, ub]
.
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.
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
— FunctionDesignMatrix(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 iteratorn
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 setout
.T
is theeltype
of the point sets. For some QMC methods (Faure, Sobol) this can beRational
It is now compatible with all scrambling methods and shifting. One can also use it with Distributions.Sampleable
or RandomSample
.
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.