Change the random number generator in Matlab function

Question

I have a task to complete that requires quasi-random numbers as input, but I notice that the Matlab function I want to use does not have an option to select any of the quasi generators I want to use (e.g. Halton, Sobol, etc.). Matlab has them as stand alone functions and not as options in the ubiquitous 'randn' and 'rng' functions. What MatLab uses is the Mersenne Twister, a pseudo generator. So for instance the copularnd uses 'randn'/'rng' which is based on pseudo random numbers....

Is there a way to incorporate them into the rand or rng functions embedded in other code (e.g.copularnd)? Any pointers would be much appreciated. Note; 'copularnd' calls 'mvnrnd' which in turn uses 'randn' and then pulls 'rng'...

• mikkola: As a first step, check out the documentation on rng.
• Matthew Gunn: This goes in the "HORRIBLE HACK, LIKELY TO CAUSE HORRIBLE PROBLEMS IN THE FUTURE" category, but if copularnd calls mvnrnd, I think you can create your own mvnrnd function in the local directory, and if you execute copularnd from that directory, your mvnrnd function will be called (instead of Matlab builtin function)...
• Allan D: Yeah @MatthewGunn. I am thinking along those lines. Note...the qrandstream and qrand/rand combo may solve the problem without causing other problems down the line. I will test it out today (once I finish other work I have to do). See my previous comment.

Answer 1

First you need to initialize the haltonset using the leap, skip, and scramble properties. You can check the documents but the easy description is as follows:

• Scramble - is used for shuffling the points
• Skip - helps to exclude a range of points from the set
• Leap - is the size of jump from the current selected point to the next one. The points in between are ignored.

Now you can built a haltonset object:

p = haltonset(2,'Skip',1e2,'Leap',1e1);
p = scramble(p,'RR2');

This makes a 2D halton number set by skipping the first 100 numbers and leaping over 10 numbers. The scramble method is 'PR2' which is applied in the second line. You can see that many points are generated:

p = 
Halton point set in 2 dimensions (818836295885536 points)
Properties:
 Skip : 100
 Leap : 10
 ScrambleMethod : RR2

When you have your haltonset object, p, you can access the values by just selecting them:

x = p(1:10,:)

Notice: So, you need to create the object first and then use the generated points. To get different results, you can play with Leap and Scramble properties of the function. Another thing you can do is to use a uniform distribution such as randi to select numbers each time from the generated points. That makes sure that you are accessing uniformly random parts of the dataset each time.

For instance, you can generate a random index vector (4 points in this example). And then use those to select points from the halton points.

>> idx = randi(size(p,1),1,4)
idx =
 1.0e+14 *
 3.1243 6.2683 6.5114 1.5302
>> p(idx,:)
ans =
 0.5723 0.2129
 0.8918 0.6338
 0.9650 0.1549
 0.8020 0.3532

• Allan D: Thanks for your input folks, but I don't think my question is adequately answered. To be more specific, I am trying to use the copularnd function in MatLab. However, it uses rng/randn (a pseudo rng, Mersenne Twister) in its code to compute the copula distribution function. I don't want to do that, I want to use quasi random numbers (Halton, Sobol etc). The documentation does not explicitly show that this is possible, thus my question. I can generate the Halton sequence...but how do I get the "rand" to use it in copularnd? @mikkola
• NKN: Then you need to explain more by editing your question and consider adding some code, or you can ask another question.
• Allan D: Update: It seems I may have to recode the copularnd function...finding a way for it to use Halton #s. For example, I may have to replace "rand" with "p" anywhere I find it in the copularnd code...I was hoping for a more elegant and efficient solution, with minimal chance of coding error.
• Allan D: The code for copularnd is as follows (uses mvnrnd) case 'gaussian' Rho = varargin{1}; n = varargin{2}; d = size(Rho,1); if isscalar(Rho) if ~(-1 < Rho && Rho < 1) error(message('stats:copularnd:BadScalarCorrelation')); end Rho = [1 Rho; Rho 1]; d = 2; elseif any(diag(Rho) ~= 1) error(message('stats:copularnd:BadCorrelationMatrix')); end % MVNRND will check that Rho is square, symmetric, and positive semi-definite. u = normcdf(mvnrnd(zeros(1,d),Rho,n))
• Allan D: That for mvnrnd is [link] (mathworks.com/help/stats/mvnrnd.html) ...etc for all the other dependancies

Answer 2

link

'qrandstream' may be the answer I am looking for....with 'qrand' instead of 'rand'

e.g..from MatLab doc

p = haltonset(1,'Skip',1e3,'Leap',1e2);
p = scramble(p,'RR2');
q = qrandstream(p);
nTests = 1e5;
sampSize = 50;
PVALS = zeros(nTests,1); 
for test = 1:nTests
 X = qrand(q,sampSize);
 [h,pval] = kstest(X,[X,X]);
 PVALS(test) = pval;
end

I will post my solution once I am done :)