RandomizedQuasiMonteCarlo

This package is now deprecated. It has been "merged" with the QuasiMonteCarlo.jl package with PR57 All the randomization methods are now available there see the doc.

The purpose of this package is to provide randomization method of low discrepancy sequences.

So far only nested uniform scrambling, Cranley Patterson Rotation (shift mod 1) and Linear Matrix Scrambling.

Compared to over Quasi Monte Carlo package the focus here is not to generate low discrepancy sequences (ξ₁, ..., ξₙ) (Sobol', lattice, ...) but on randomization of these sequences (ξ₁, ..., ξₙ) → (x₁, ..., xₙ). The purpose is to obtain many independent realizations of (x₁, ..., xₙ) by using the functions shift!, scrambling!, etc. The original sequences can be obtained for example via the QuasiMonteCarlo.jl package.

The scrambling codes are inspired from Owen's R implementation that can be found here.

Basic examples

```julia using RandomizedQuasiMonteCarlo, QuasiMonteCarlo m = 7 N = 2^m # Number of points d = 2 # dimension M = 32 # Number of bit to represent a digit b = 2 # Base uuniform = rand(N, d) # i.i.d. uniform usobol = permutedims(QuasiMonteCarlo.sample(N, zeros(d), ones(d), SobolSample())) # I should update the convention in my pkg to have dim × n and not n × dim unus = nesteduniformscramble(usobol; M=M) ulms = linearmatrixscramble(usobol, b; M=M) udigitalshift = digitalshift(usobol, b; M=M) ushift = shift(usobol)

Plot

using Plots, LaTeXStrings

Settings I like for plotting

default(fontfamily="Computer Modern", linewidth=1, label=nothing, grid=true, framestyle=:default)

begin d1 = 1 d2 = 2 sequences = [uuniform, usobol, unus, ulms, ushift, udigitalshift] names = ["Uniform", "Sobol (unrandomized)", "Nested Uniform Scrambling", "Linear Matrix Scrambling", "Shift", "Digital Shift"] p = [plot(thicknessscaling=2, aspect_ratio=:equal) for i in sequences] for (i, x) in enumerate(sequences) scatter!(p[i], x[:, d1], x[:, d2], ms=0.9, c=1, grid=false) title!(names[i]) xlims!(p[i], (0, 1)) ylims!(p[i], (0, 1)) yticks!(p[i], [0, 1]) xticks!(p[i], [0, 1]) hline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2) vline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2) hline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8) vline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8) end plot(p..., size=(1500, 900)) end ```

Scrambling of Faure sequences in base b

Now let say you want to do scrambling in base b.

julia using QuasiMonteCarlo using Random

The InertSampler is used to obtain a randomized Faure sequences (see here in the test file).

julia struct InertSampler <: Random.AbstractRNG end InertSampler(args...; kwargs...) = InertSampler() Random.rand(::InertSampler, ::Type{T}) where {T} = zero(T) Random.shuffle!(::InertSampler, arg::AbstractArray) = arg

```julia m = 4 d = 3 b = QuasiMonteCarlo.nextprime(d) N = b^m # Number of points M = m rng = InertSampler()

Unrandomized low discrepency sequence

u_faure = permutedims(QuasiMonteCarlo.sample(N, d, FaureSample(rng)))

Randomized version

unus = nesteduniformscramble(ufaure, b; M=M) ulms = linearmatrixscramble(ufaure, b; M=M) udigitalshift = digitalshift(ufaure, b; M=M)

```

This plot checks (visually) that you are dealing with $(t,d,m)$ sequence i.e. you must see one point per rectangle.

```julia begin d1 = 1 d2 = 3 x = ulms p = [plot(thicknessscaling=2, aspect_ratio=:equal) for i in 0:m] for i in 0:m j = m - i xᵢ = range(0, 1, step=1 / b^(i)) xⱼ = range(0, 1, step=1 / b^(j)) scatter!(p[i+1], x[:, d1], x[:, d2], ms=2, c=1, grid=false) xlims!(p[i+1], (0, 1.01)) ylims!(p[i+1], (0, 1.01)) yticks!(p[i+1], [0, 1]) xticks!(p[i+1], [0, 1]) hline!(p[i+1], xᵢ, c=:gray, alpha=0.2) vline!(p[i+1], xⱼ, c=:gray, alpha=0.2) end plot(p..., size=(1500, 900)) end

```

Multiple randomization

In case you need to repeat randomization several times, I suggest you use in place functions and compute some stuff in advance e.g. bit expansion of your initial set of point to be randomized, the which_permutation function used in Nested Uniform Scrambling.

julia unrandomized_bits = sobol_pts2bits(m, d, M)

Here I use directly the bit representation for Sobol sequence. I could have done like in previous example and import the Sobol sequence from QuasiMonteCarlo.jl (which calls Sobol.jl). Note that you can call any sequence of point into its bit representation with points2bits function.

```julia randombits = similar(unrandomizedbits) indices = whichpermutation(unrandomizedbits) # This function is used in Nested Uniform Scramble. I nus = NestedUniformScrambler(unrandomizedbits, indices) lms = LinearMatrixScrambler(unrandomizedbits)

usob = dropdims(mapslices(bits2unif, unrandomizedbits, dims=3), dims=3) unus = copy(usob) ulms = copy(usob) ushift = copy(usob)

NumberOfRand = 100

for j in 1:NumberOfRand scramble!(unus, randombits, nus) scramble!(ulms, randombits, lms) shift!(u_shift) end

```