Quasi-MH

I think that taking advantage of this should let you get a significant improvement in sampling without being too much work.

[just-cameron]

Oh, interesting! Hadn't seen that before. Might be worth implementing when I get a moment.

[devmotion]

Is there actually anything to implement for AdvancedMH here? I skimmed through the paper and it seems that the user just has to provide an RNG that returns a QMC sequence. There are established packages that generate QMC sequences, so I don't think that they should be reimplemented in AdvancedMH. And it is already possible to use custom RNGs in AdvancedMH (and generally any AbstractMCMC dependent package).

Is there actually anything to implement for AdvancedMH here? I skimmed through the paper and it seems that the user just has to provide an RNG that returns a QMC sequence. There are established packages that generate QMC sequences, so I don't think that they should be reimplemented in AdvancedMH. And it is already possible to use custom RNGs in AdvancedMH (and generally any AbstractMCMC dependent package).

I'm uncertain about the internal implementation details of AdvancedMH; it depends on how AdvancedMH uses random numbers. For instance, let's say we get a user-provided van der Corput generator seeded to begin at .25. This returns 1/4, 3/4, 1/8, 3/8, 5/8, 7/8... Say the sampler first draws one value to determine the acceptance threshold, then draws one value to propose a new point. This will give inconsistent results -- note that there is a strong correlation between n and n+1, which tend to be similarly sized. Thus, proposed points below 1/2 will be rejected often, and proposed points greater than 1/2 will hardly be rejected at all.

Instead, you'd need to generate some kind of multidimensional low-discrepancy sequence which is completely uniformly distributed, like a generalized golden ratio sequence, which guarantees very low correlations between the proposals and acceptance probabilities (lower than you would get by chance, even).

I'm also unsure how computation of e.g. standard errors or MCSE will work; this might require a different calculation, since QMC is better than normal MCMC by a factor of sqrt(n).

[devmotion]

Sure, this example is also mentioned in the paper. As far as I understood, therefore the authors propose to use so-called CUD sequences.

In general, it seems difficult to avoid this problem if users provide a "bad" RNG since it is the number of (quasi-)random samples that are used to generate the proposal is not fixed (e.g., generally a proposal from a multivariate normal distribution of dimension d uses d (quasi-)random samples).

@devmotion The main issue I'm referring to is that the implementation of a quasi-MH needs to take in a random number generator that returns an array of uniformly-distributed random numbers (one for the acceptance threshold, plus one for each dimension), rather than a single random integer like a Mersenne Twister gives. Drawing every number from a single low-dispersion sequence, rather than a low-dispersion sequence in multiple dimensions (Like a Sobol sequence) will generate correlations between proposal points and acceptance probabilities, when the goal of quasi-MH is to minimize these correlations.

[devmotion]

As I said, I am aware of the problems. My main point is that it seems difficult to ensure that these issues are avoided in the abstract interface in AdvancedMH since we don't have enough information about the proposal. However, the user does, and hence my suggestion to make it the user's responsibility to avoid these problems.

As I said, I am aware of the problems. My main point is that it seems difficult to ensure that these issues are avoided in the abstract interface in AdvancedMH since we don't have enough information about the proposal. However, the user does, and hence my suggestion to make it the user's responsibility to avoid these problems.

I'm confused. Does AdvancedMH let you provide an RNG that returns a vector of random numbers, rather than a single random number, then?

[devmotion]

AdvancedMH let's you provide an arbitrary RNG. And in particular you could provide a custom RNG that returns quasi-random numbers in the required order and with the required properties.

AdvancedMH let's you provide an arbitrary RNG. And in particular you could provide a custom RNG that returns quasi-random numbers in the required order and with the required properties.

Alright, I think that makes some sense, but I don't think it's at all intuitive from the user's perspective how to pass something like a Sobol sequence to the RNG without prying open the code and figuring out exactly how the code uses the provided random number generator. Anything else risks mangling the structure and getting inconsistent results, if it turns out that AdvancedMH uses random numbers differently than expected. For instance, does AdvancedMH require arrays of random numbers, or just one random number at a time? If so, how do I know which random sequence will be used for which dimension, and which is used for acceptance thresholds? What if I try to generate multiple chains?

[devmotion]

For instance, does AdvancedMH require arrays of random numbers, or just one random number at a time? If so, how do I know which random sequence will be used for which dimension, and which is used for acceptance thresholds?

The whole point I want to make is that AdvancedMH doesn't know any of these things. It just passes the RNG to the user-provided proposal and uses randexp to sample an acceptance threshold. AdvancedMH does not know how many random numbers are required and it does not even sample anything explicitly apart from the acceptance threshold. To me it seems the only person who can possibly know how and in which order random numbers are generated is the user.

What if I try to generate multiple chains?

The default sampling algorithms in AbstractMCMC for multiple chains with multiple threads or processes use independent copies of the user-provided RNG (with randomly initialized reproducible seeds).

The whole point I want to make is that AdvancedMH doesn't know any of these things. It just passes the RNG to the user-provided proposal and uses randexp to sample an acceptance threshold. AdvancedMH does not know how many random numbers are required and it does not even sample anything explicitly apart from the acceptance threshold. To me it seems the only person who can possibly know how and in which order random numbers are generated is the user.

I think I kind of understand your point, but my point is that for a user, trying to figure out how to implement quasi-MC by fiddling with the RNG is going to be difficult and will very likely lead to inconsistent estimates. Another example of this is that I have no idea how the distributions I pass will end up implementing the by using the random number I give them. For instance, if normal variables are generated in pairs using the Box-Muller transform, this can also result in errors when using quasirandom sequences, since the Gaussian proposals will no longer be independent. I think it makes far more sense to just implement a single quasirandom sequence that works well out of the box, rather than asking people to comb through the source code to provide an appropriate sequence and risk many of them accidentally creating inconsistent estimators.