Created
April 20, 2021 21:33
-
-
Save cscherrer/97305d62f83b02dec62a86d4916305c5 to your computer and use it in GitHub Desktop.
`rand` method for "hypercubes", e.g. Sobol sequence
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using Sobol | |
using MeasureTheory | |
abstract type Hypercube{k} end | |
struct SobolHypercube{k} <: Hypercube{k} | |
seq :: SobolSeq{k} | |
value :: Vector{Float64} | |
index :: Ref{Int} # start at zero | |
function SobolHypercube(k::Int) | |
seq = SobolSeq(k) | |
value = Sobol.next!(seq) | |
return new{k}(seq, value, 0) | |
end | |
end | |
function Sobol.next!(s::SobolHypercube) | |
s.index[] = 0 | |
Sobol.next!(s.seq, s.value) | |
end | |
function rand(ω::SobolHypercube{k}) where {k} | |
ω.value[ω.index[] += 1] | |
end | |
function rand(ω::Hypercube, d::ParameterizedMeasure) | |
Dists.quantile(distproxy(d), rand(ω)) | |
end | |
function rand(ω::Hypercube, d::ProductMeasure) | |
[rand(ω, dj) for dj in d.data] | |
end | |
using Colors | |
using GLMakie | |
function makeplot() | |
ω = SobolHypercube(100) | |
x = range(-2,2,length=100) | |
d = For(x) do xj Normal(2*xj, 1/(1 + xj^2)) end | |
y = rand(ω, d) | |
plt = scatter(x,y,markersize=2, color=colorant"rgba(0,0,0,0.1)") | |
for j in 1:100 | |
next!(ω) | |
y = rand(ω, d) | |
scatter!(x, y, markersize=2, color=colorant"rgba(0,0,0,0.1)") | |
end | |
plt | |
end | |
makeplot() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment