Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sobol' Sequence Missing First 4 Points #78

Closed
alegresor opened this issue May 17, 2023 · 2 comments
Closed

Sobol' Sequence Missing First 4 Points #78

alegresor opened this issue May 17, 2023 · 2 comments

Comments

@alegresor
Copy link

Running the following code

using QuasiMonteCarlo
QuasiMonteCarlo.sample(4,0,1,SobolSample())

yields

4-element Vector{Float64}:
 0.375
 0.875
 0.625
 0.125

However, the first 4 points in a 1 dimensional Sobol' sample should be [0, .5, .25, .75]. The output are points at index 5 to 8.
I am looping in my collaborator @fjhickernell, an expert in QMC methods who will be valuable to this discussion.

@dmetivie
Copy link
Contributor

I believe this is caused by the “famous” debate about starting the Sobol sequence at zero (-> theory) or not (-> practical implementation).
QuasiMonteCarlo.jl uses Sobol.jl which use the latter convention (drop zero), see related issue and refs therein (Owen' paper).
This debate was started way back I believed and amplified with Pytorch using the drop zero solution. @fjhickernell probably knows better.

Meanwhile, a while ago, I ported the R Owen (@dcarrera you suggested that at some point) code to do scrambling in a package RandomizedQuasiMonteCarlo.jl.
There, you can find the implementation starting at zero.

Note: I don't maintain the package anymore since all its features have been moved in QuasiMonteCarlo.jl (excepted the following function, in fact).

using RandomizedQuasiMonteCarlo
m = 2
N = 2^m # Number of points
d = 1 # dimension
b = 2 # always 2 for Sobol
M = 32 # number of bits to encode points
unrandomized_bits = sobol_pts2bits(m, d, M) 
u_sobol = dropdims(mapslices(bits2unif, unrandomized_bits, dims=3), dims=3)

yields

4×1 Matrix{Rational{Int64}}:
 0//1
 1//2
 1//4
 3//4

I believe this convention should be specified somewhere in bold in Readme.md and docs.

@alegresor
Copy link
Author

Thank you for the quick response. I'll close this issue and monitor the other one you linked.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants