Quasi Monte Carlo in Matlab

Question

I want to use Quasi Monte Carlo to try and improve the convergence of a simulation I am running.

The random numbers are simply to produce the observation errors for a standard linear regression model, which is then estimated using a number of different regression techniques. This is done repeatedly to estimate the mean square error of each model.

I'm fairly new to Quasi Monte Carlo, but is is likely to help in this situation I am just using it to produce 10k random numbers. It seems that generally I can expect quicker convergence of the order of $1/n$ rather than $n^{-1/2}$: http://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method.

However, the above article also states that the QMC numbers are not truly random, so I wonder what the implications might be for any statistical tests I might want to run on the results.

1.) What are the pros and cons of MC vs. QMC. (Would you always want to use QMC if it's available?) 2.) What tests can I use to ascertain which is best for my application? (Seems any test that depends on the numbers being truly random will fail?)

I know that this can be done in Matlab using

q = qrandstream('halton',NSteps,'Skip',1e3,'Leap',1e2); 
RandMat = qrand(q,NRepl); 
z_RandMat = norminv(RandMat,0,1);

which is taken from this paper.

There are other low discrepancy numbers, such as Sobol sequence also available in Matlab, and again I'd like to know what tests I can use to ascertain which is best for my situation.

Answer 1

Quasi random numbers are not statistically independent so if your algorithm requires this, you should not use them.

One area where you can usually use quasi random numbers instead of pseudo random numbers is in the evaluation of integrals via monte carlo integration. In this application, the fact that quasi random numbers are not truly random is not an issue since the statistical properties don't matter. What matters is how well the sequence samples the integrand space and quasi random numbers fill space better than pseudo random numbers, leading to improved convergence.

I find the following article by the Numerical Algorithms Group (NAG) useful in describing this in more detail

http://www.nag.co.uk/IndustryArticles/introduction_to_quasi_random_numbers.pdf

I also note that Mathworks suggest a couple of other applications of quasi-random numbers along with monte-carlo integration: space filling experimental designs and global optimisation http://www.mathworks.co.uk/help/stats/generating-quasi-random-numbers.html#br5k9hi-9

One issue you need to be wary of is how many dimensions a quasi-random sequence retains its desirable properties for. The standard Halton sequence, for example, does not do well in high dimensions. You can see this directly using MATLAB

Create a 15 dimensional halton sequence and evaluate 1000 points in each dimension:

>> rng default
>> p=haltonset(15);
>> x0=net(p,1000);
>> size(x)

ans =
        1000          15

Plot dimensions 1 and 2 against each other

>> scatter(x0(:,1),x0(:,2),5,'b')

Looks good so far. Not so good when we plot dimension 14 against 15

>> scatter(x0(:,14),x0(:,15),5,'b')

There are various ways of dealing with this such as scrambling (see http://www.mathworks.co.uk/help/stats/qrandset.scramble.html for example) but you may also choose to select different sequences that are known to be good for higher dimensions.

I have more experience with using the NAG Toolbox for MATLAB (http://www.nag.co.uk/numeric/MB/start.asp) which provides 2 sobol generators (one that's good for 10,000 dimensions, another thats good for 1111 dimensions), a Niederreiter generator (318 dimensions) and a Faure generator (40 dimensions).