How to determine how many random numbers were generated, using Mersenne twister?

I want to use Mersenne Twister to generate 2 sets of random numbers using common seed, and want to determine which of these sets used more random numbers. Random number generation is performed by a called function; so calling function doesn't directly know how many random numbers were generated in each call.

rng_state_initial = rng;
% Generate 1st sequence
[output1] = joeblow1(...);
% save the state of rng after 1st sequence generated
rng_state_after_1st_sequence = rng;
% Generate 2nd sequence
rng(rng_state_initial);
[output2] = joeblow2(...);
% save the state of rng after 2nd sequence generated
rng_state_after_2nd_sequence = rng;

If I want to then generate additonl lrandom numbers indepedent of 1st and 2nd sequences, should I use

rng(rng_state_after_1st_sequence)

or

rng(rng_state_after_2nd_sequence)

before generating the new random numbers?

I don't know which random number state (1st or 2nd sequence) is further along.

I can do a hack job to do either rng(rng_state_after_1st_sequence) or rng(rng_state_after_2nd_sequence) and then generate a ginormous number of random numbers that I'm sure exceeds the number of random numbers generated in 1st or 2nd sequence,, save the state after that, and then generate the new random numbers. But I don't want to do that.

Note that multiple streams are not available with Mersenne Twister.

[Jeff Miller]

Another cluges, if you know an upper limit on the difference between the numbers of rng calls by the two functions. E.g., the one that is farther ahead is at most 1000 calls ahead of the one that is farther behind. If so, it seems like you could generate 1000 numbers from sequence 1 and see whether it catches up to the sequence 2 (i.e., generates the next number produced by sequence 2). You could then proceed using whichever sequence you determine to be farther ahead.

If you know an upper bound like 1000 that isn't too large, it might be a bit faster than the cluge you are thinking of.

[Mark Stone]

Yes, my hack job accounts for stuff like this. That will work on my existing function. It could break down on some other user's new function, which might be way out in left field vs. what I anticipated; hence why the hack job is unsatisfying for use by someone other than me. If I don't get a "nice" answer, I'll implement a hack job for now, while I consider a bettter long term solution.

Answer 1

[Walter Roberson]

Record the two rng states. Loop generating one number at a time from each of the two streams until the newly generated state of one of them equals the recorded state of the other.

There is no iteration count information recorded in the state.

If I recall correctly, about two years ago I saw an abstract for a paper about recovering rng state. If you generate N random numbers, can you deduce what the rng state must have been? The number required was much smaller than you might expect. I do not recall the details now, but my brain is saying "Wasn't it only in the range of 8 to 13 needed?" I did not read the paper to see what techniques they were using, but I wonder if some of what they did might be useful for figuring out the meaning of the state differences?

[Mark Stone]

Thanks. Unfortunately, this seems like it will be painful in run time. But I'll give it a go.

Answer 2

[Jan]

The direct solution is to modify the function joeblow1 and joeblow1 such, that they count the created random numbers.

Another option is to use a wrapper function, which calls the random number generation and counts the number of created values.

Using the state of the random number generator to count the produced values is rather indirect and cannot be efficient. I cannot imagine a purpose to create additional random numbers only to set the rng to a defined state. Using a new state would be much cheaper. Therefore I assume, that this is an XY-problem. Please explain, which problem you want to solve actually.

"Note that multiple streams are not available with Mersenne Twister." - is this really the case?

[Mark Stone]

"Note that multiple streams are not available with Mersenne Twister." - is this really the case?

Yes,per the below:

https://www.mathworks.com/help/matlab/math/creating-and-controlling-a-random-number-stream.html

Choosing a Random Number Generator

mt19937arMersenne twister entry for Multiple Stream and Substream Support is "No"

++++++++++++++++++++++++++++++++++++++

I have a large library of joeblow functions which are called by multiple functions, and wish to retain backward compatibility, and not mess with changing interfaces. When all this was previously devleoped, I didn't see the need to count random numbers. Now I am trying to add a new mode of operattion to my program which needs common random numbers in each "iteration", and independent random numbers across iterations. And for various reasons, I wish to keep using Mesenne Twister.

If no solutiion is offered, I will resport to a hack job. My hack job woudl generate a bunch of new random numbers, which would be "thrown away", and only used to advance the random number generator state to where subsequent random numbers would be independent of the previous.

I was hoping there is a way of determining which of two random number gnerastor states is further along (or how many random numbers were generated vs. previous random number genrator state).

[Jan]

I'd solve this with a wrapper:

function out = myRand(varargin)
persistent count
if nargin == 1 && strcmpi(varargin{1}, 'flush')
    out   = count;
    count = 0;
    return;
end
   
out = rand(varargin{:});
end

Then all calls to rand() inside blowjoe have to be changed in the editor.

I'm using my own random number generators written in C, which allow to include a counter easily.

[Mark Stone]

Unfortunately, this violates my goal of backward compatibility. I'd rather do the hack job than have such a solution. If I were just using the program myself, the hack job would work. But not so pretty or bulletproof if I ever release it for other peope to use.

I'm still hoping somone knows how to extract what I need from the state returned by rng;

[Jan]

You can simply include myRand() in the code you distribute. This is even backward compatible to the Matlab versions before rng was introduced.