rng state appears to be dependent on distribution type

I'm converting some random number generation stuff from python to matlab and I'm noticing some strange behavior in how the rng updates its state. It appears that in python the random number generator state is only dependent on the number of random draws, but in matlab its dependent on the number of draws and the distribution type. The consequence of this is if you are generating a sequence of random draws (ie a Monte Carlo) and you change the distribution type for one item in the middle of the sequence it ends up changing the random numbers you get for all the random draws after that point. Does anyone know if this behavior is correct? I'm used to things behaving more like how python is showing.

Python check:

import random 
#% Perform 1000 uniform draws
random.seed(1) 
draw1 = [random.uniform(0,1) for x in range(1000)] 
state1 = random.getstate()
#% Repeat same sequence, but with 1000 gaussian draws
random.seed(1) 
draw2 = [random.gauss(0,1) for x in range(1000)] 
state2 = random.getstate() 
#% Check if rng ends at same state for both tests
print(state1 == state2) 
>> True

Versus matlab:

% Perform 1000 uniform draws
rng(1) 
draw1 = rand(1000,1); 
state1 = rng; 
% Repeat same sequence, but with 1000 gaussian draws
rng(1) 
draw2 = randn(1000,1);
state2 = rng; 
% Check if rng ends at same state for both tests
disp(isequal(state1,state2)) 
>> 0

[Paul]

I think it depends on the algorithm used by randn. For example, some algorithms use more than one random number, e.g., from a uniform distribution, to generate a standard normal random number.

[Paul]

Sometimes, changing the parameter of the distribution can change the sequence of random draws, at least for mvnrnd, as discussed in this thread.

Answer 1

[Walter Roberson]

randn() has different ways of computing a normal sample, with the exact method depending on which variety of random number generator is being used. The default random number generators use a ziggurat method of sampling for randn()

https://blogs.mathworks.com/cleve/2015/05/18/the-ziggurat-random-normal-generator/

I have limited internet access at the moment so I cannot presently find the reference; I seem to recall that the algorithm on average consumes about 1.3 "rand" equivalent

[Steven Lord]

"About 1.3 uniform values" is if you specify the 'Polar' NormalTransform method on a RandStream object. If you use 'Inversion' that requires exactly 1 uniform value. 'Ziggurat' (the default) requires about 2.02 uniform values per normal value. See the "Transformation Algorithms" section on this documentation page.

rng(1);
rand(1, 1000);
s1 = rng;

rng(1);
gs = RandStream.getGlobalStream;
gs.NormalTransform = 'Inversion';
randn(1, 1000);
s2 = rng;

isequal(s1, s2)

ans = logical
   1

[Michael]

Thanks was unaware that there were different underlying algorithms for a simple normal distribution.

[Steven Lord]

Sometimes things that seem simple have many, many devils when you look closer at the details :)