Implementation of scrambled Van der Corput sequence differs from the referenced paper

[treverhines]

I noticed a potential issue with the function scipy.stats._qmc.van_der_corput. That function is based on Algorithm 1 in this paper:

https://arxiv.org/pdf/1706.02808.pdf

However, there is a difference between the scipy implementation and the algorithm from that paper. In Algorithm 1, the last line in the while loop is res <- (res - dig) / b, but the scipy implementation is effectively doing res <- (res - pdig) / b, using the variable names from Algorithm 1. This difference does amount to different output sequences.

Since dig is res mod b, the line res <- (res - dig) /b is just integer division. The scipy implementation can then be made consistent with the referenced algorithm by changing this line:

scipy/scipy/stats/_qmc.py

Line 522 in f856e3e

quotient = (quotient - remainder) / base

to

quotient //= base

I am currently working on cythonizing van_der_corput (see the discussion here: #13474), so I can also make the change in that PR (which I will probably make today).

[tupui]

Thanks for reporting. But not sure this would work for the scrambled case as there is a permutation happening on the remainder.

FYI, a good way to check that you did not break anything is to look at the convergence plots in stats' tutorial.

[treverhines]

Sorry for not being more clear. I am saying that the algorithm in Owen 2017 (which is for the scrambled case) is subtracting the unpermuted remainder in the last line of the while loop, but the current implementation is subtracting the permuted remainder.

Thanks for pointing out the convergence plots. I will check that out.

[tupui]

Oh ok no my bad I did not click the reference you pointed out. I did this with Art. I will dig up my emails and see if we talked about that

[tupui]

@treverhines just checking in 😃 Did you get some time here?

Looking a bit more at what you said. I think you are right even if it seems to not have a huge impact. I have quickly redone the convergence plots following what I did back then (https://gist.github.com/tupui/ee6e685219c06c429338569e7ee7ff2a) and we have the right convergence rate. Although I would need to get both versions on the same curve to be sure. Anyway, thanks again for the finding.

[tupui]

I did the plots and there is basically no change. So feel free to do it in the same PR as it would be closer to the published algorithm 😃

Just showing one the results using my script here https://gist.github.com/tupui/ee6e685219c06c429338569e7ee7ff2a